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ABSTRACT 


The Navier-Stokes equations are solved by numerical 
method for steady, fully developed, incompressible, laminar 
flow in curved rectangular channels considering the curvature 
ratio effect in the formulation. Solutions are obtained for 
aspect ratios 1, 2, 5 and 0.5 and Dean number ranges from 5 
to 715, for example, for the case of curved square channel. 

It is found that an additional counter-rotating pair of vortices 
appears near the central outer region of the curved channel in 
addition to the familiar secondary flow at a certain higher Dean 
number. The correlation equations for the friction factor are 
developed. The friction factor results are compared with the 
available theoretical and experimental results and the agreement 
is found to be good. The curvature ratio effect is studied and 
it is confirmed that when r > 10, Dean's original formulation 


neglecting curvature ratio effect is applicable. 


The Graetz problem in curved square channels considering 
the curvature ratio effect in the formulation is studied for the 
two basic thermal boundary conditions of constant wall temperature 
and uniform wall heat flux using the accurate velocity fields. 

The roles of Prandtl] and Dean numbers on heat transfer results 


are found to be similar. When Prandtl or Dean number is large, 
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the local Nusselt number fluctuates before the asymptotic value 
is reached. The asymptotic Nusselt numbers are of practical 
importance and a simple correlation equation valid for the 
cases of constant wall temperature and uniform wall heat flux 


for curved square channels is developed. 


The Graetz problem in curved rectangular channels 
with uniform convective boundary condition is also studied. 
The cases of Biot numbers 1 and 10 are studied in detail. 
The secondary flow effect on liquid solidification-free zone in 
curved rectangular channels is investigated for Pr = 0.03 (liquid 


metal) and Pr = 10.05 (water). 
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CHAPTER I 


INTRODUCTION 


Curved channels for fluid flow are used in various 
engineering applications, notably in heat exchanger and air 
conditioning apparatus, In contrast to fluid flow in the 
Straight channels, curved channel flows are characterized by 
secondary flow in the cross section normal to the mainflow 
as a result of centrifugal forces acting on the fluid. 

The phenomenon of steady fully developed laminar 
flow of a viscous fluid in a curved tube is now well understood 
and accurate friction factor correlation equations are available 
in the literature. For the analytical solution of the problem, 
Dean's perturbation analysis [1,2] is valid for small Dean 
number flow regime and the asymptotic boundary-layer analysis 
[3,4] is applicable only for large Dean numbers. Recently, 
numerical solution [5,6] has been carried out for Dean number, 
Re(a/R) !/2, ranging from 1 to 1200 with curvature ratio parameter 


(R/a) varying from 10 to 100. 


It is expected that the method of the solution for 
the equally important case of fully developed laminar flow in 
curved rectangular channels is similar to that for curved tubes. 


Indeed, the perturbation solutions [7,8] valid for small Dean 
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numbers and asymptotic solutions [9,10] for curved square 
channels .based on boundary-layer approximation and applicable 
to large Dean numbers are available. The flow results from 
numerical solutions [11] are reported for curved rectangular 
channels with aspect ratios y = 0.2, 0.5, 1, 2 and 5 up to Dean 


1/2 


number, K = Re(De/R) ‘“, of order 100. 


Heretofore, the theoretical analysis for flow in 
curved rectangular channels is limited to Dean's original 
formulation where the curvature ratio R/De is assumed to be 
large and the curvature ratio parameter does not appear 
independently in the governing equations. In view of the approx- 
imate nature of the boundary-layer method and the limited 
applicability of the perturbation solution, it is desirable 
to obtain numerical solution for Navier-Stokes equations including 
the curvature-ratio effect which govern the steady fully developed 


laminar flow in curved rectangular channels. 


The superposition of the secondary flow on the axial 
flow has the effect of orderly mixing, which is expected to 
increase the heat and mass transfer rates. An effort will be 
made to bring out the heat transfer mechanism in the thermal 
entrance region. Actually, laminar forced convection heat transfer 
in the thermal entrance of curved pipes or channels is an extension 


of the classical Graetz problem and is one of the basic problems 
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in the convective heat transfer. 


It is only recently that the Graetz problem has been 
studied for the case of curved circular pipes [16,19,20,21]. 
Heat transfer results in curved rectangular channels are 
reported so far only for the limiting conditions corresponding 
to the hydrodynamically and thermally fully developed flow fields 
[11,14]. With the assumption of constant physical properties, 
the momentum and energy equations are uncoupled and the Navier- 
Stokes equations governing a steady fully developed laminar flow 


can be solved independently. 


iy Slag ts seetab ite to study the Graetz problem in curved 
rectangular channels based on accurate flow field obtained from 
the numerical solution of Navier-Stokes equations considering 
the curvature ratio effect in the formulation. In order to 
limit the scope of present study, the Graetz problem in curved 
Square channels is studied for the two basic thermal boundary 
conditions of constant wall temperature and uniform wall heat 


flux and uniform convective boundary condition. 


In contrast to the first two boundary conditions, the 
available solutions [25, 26, 27, 28, 29, 30] to Graetz problem 
with convective boundary condition are fairly limited and are 
so far confined to straight tubes and parallel-plate channel 


only. No investigation has been reported for Graetz problem in 
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curved pipes or channels with convective boundary condition in 


spite of its practical and theoretical importance. 


It should be pointed out that with convective boundary 
condition the solidification of liquids flowing in the curved 
tubes or channels may occur under extreme ambient conditions. 

In practice, the external heat transfer coefficient may vary 
circumferential ly or axially. However, for simplicity, the 
external convection will be assumed to be uniform cooling. 

The present study based on uniform convective cooling yields 

the length of solidification-free zone preceding the region 

where the solidification Occurs. The cases of water and liquid 
metal are studied in detail because of their practical importance. 
The aspect ratio effect is studied in the solidification-free 


zone analysis also. 


In summary, the objectives of this investigation are: 
1. Presentation of numerical results for fully developed laminar 
flows in curved rectangular channels with aspect ratios y = 0.5, 
1, 2 and 5 up to higher Dean number flow regime. 
2. The numerical solution of Graetz problem in curved square 
channels for the two basic thermal boundary conditions of constant 
wall temperature and uniform wall heat flux. 
3. The numerical solution of Graetz problem in curved rectangular 


channels for uniform convective boundary condition and a study of 
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the effect of secondary flow on solidification-free zone for liquid 


metal and water. 


An effort will be made to compare the numerical results 
from the present analysis with the theoretical and experimental 


results available in the literature. 
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CHAPTER II 


FULLY DEVELOPED LAMINAR FLOW IN CURVED 
RECTANGULAR CHANNELS 


2.1 Governing Equations 

For the theoretical analysis of the flow in curved 
channels,a toroidal coordinate system as illustrated in Fig. 1 is 
introduced. One assumes that the flow is steady, fully developed, 
laminar and incompressible. The detailed derivation of the 
basic governing equations is given in Appendix A. The results 
are: 
Continuity equation 
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Momentum equation in X direction 
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Momentum equation in Y direction 
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The corresponding boundary conditions are 


U=V=W=0 at walls 
= 0 along the line of symmetry Y=0 (5) 
It is convenient to introduce the following dimensionless variables: 


V= 
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The governing equations then become; 
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One can eliminate the pressure terms from the secondary flow 
equations and obtain the following vorticity transport equation 


for secondary flow. 
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The dimensionless stream function is defined as 
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Then one obtains the following stream function-vorticity equation. 
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The axial flow equation remains the same. It is restated as; 
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The boundary conditions for the above three variables are; 


(a) the axial velocity w 


w = 0 at walls 
CAR 0 along the line of symmetry y=0 (14) 


(b) the stream function y 


vw = 0 at walls and along the line of symmetry y=0 (15) 
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(c) the vorticity & (16) 


— is not known in advance at the walls. This will be 
discussed in detail in the next section, However, & = 0 along 


the line of symmetry y=0. 


At this point, it is noted that a numerical solution 
is the only practical approach for the solution of the above 
set of simultaneous quasilinear, second-order partial differential 


equations of elliptic type. 


Before the numerical solution technique is introduced, 
two alternative expressions for the product of friction factor 


and Reynolds number, fRe, are presented here: 
1. The expression based on the velocity gradient along the channel walls. 


The surface area of the upper half of a rectangular 


channel along the axial length AZ = RAd is 
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Using the non-dimensional quantities, one has 
(4) (1+y) 
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The product of friction factor and Reynold's number is defined as; 
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13 
when 1/r <<1, a simplified expression for (fRe). becomes 


s, ow = 
(FRe). = 2( Fi )//" (20a) 
II. The expression based on the overall force balance for the 
differential axial length dZ. 


A force balance on a flow element with axial length dZ 


gives; 
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These two alternative methods of calculation for (fRe) 
provide a means of assessing the accuracy of the numerical solution. 


This will be discussed in the finite difference solution section. 
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2.2 Finite Difference Solution 

Using the three-point central difference approximation 
for all the derivatives in equations (10), (12) and (13), each 
equation can be transformed into a set of linear algebraic 


equations in the following form: 
A ae sop ita pean 22k Oe ae 


where F is a dummy variable which stands for €, wy and w,and A, B, 
C and D are the coefficients with j=2,...N. For the vorticity 
transport equation (10), F = & and the coefficients A, B, C and 
D are given as: 
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For the stream function vorticity equation (12), one has Fey 


and the coefficients are 
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Similarly, with Few for axial velocity equation (13), the 


coefficients become 
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The secondary velocity components u and v are computed from the 


following finite-difference equations. 
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The boundary conditions can be written in finite difference 


form as; 
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As far as the vorticity at the boundaries is concerned, 
it is treated in the following way. One can represent the 
vorticity at the inner wall as E14? j=l,...,N+l, referring to the 
grid point system as shown in Fig. 1. At the wall, u=v=0, and 


ou/ay = 0 but av/ax # 0, So one can obtain 
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By using Taylor's series expansion, one has 
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here 1 stands for the position at wall, 2 for the interior point 
next to the wall, and 3 for the interior point next to 2. From 


the boundary conditions, one knows that 
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2 2 3 3 
vy = (DY) APD, Ay, L2OxD + ofan) (32b) 
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Neglecting the second term of equation (32a) one obtains 


(33a) 


and from (32a) and (32b), eliminating the third order term, one 


has 
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Equation (33a) has larger truncation error than equation (33b), 
but it is a stable representation. Equation (33a) is employed 
in this study. The vorticity at walls is expressed in the 
following way: 
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For the upper wall 
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Using the boundary conditions stated, the set of 
finite difference equations can be written in the following 
matrix form after applying equation (22) to the grid points along 
the line j=2,...,N 
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It is noted that a tridiagonal matrix results for the 
system of algebraic equations for each line j=2,...,N. The 
finite difference solution of the coupled equations together 


with the associated boundary conditions are carried out in the 
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following manner: 

ae For given curvature ratio r and axial pressure gradient, 
the initial values for dependent variables &, w and w are 
provided. 

ae Using the initial values for stream function y, 

the boundary vorticity is computed using equation (34, 35, and 
36). After obtaining the boundary vorticities, the vorticity 
transport equation (23) is solved by the Gaussian elimination 
method [12] along the lines j=2,...,N. 

oe The elliptic equation for stream function, wv, is solved 
by a line iterative method [12] using a relaxation factor w as 
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The iteration is terminated when the following relative error 


criterion is satisfied 
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equation (26). Using the updated secondary velocity components, 
the axial momemtum equation can be solved following exactly 
the method described in step (3). Thus one completes one cycle 


of outer iteration. 
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Sp Steps 2 to 4 are repeated until the following convergence 
criterion is satisfied for wy and w. 


(k#1) — (k) (K+) : 
max le ie | / max ieee Peevto 


3 
where k stands for kth outer iteration. 

Numerical experiments are conducted to establish the 
mesh size and the values of the relaxation factor w, ranging 
from 1 to 0.1, are used as Dean number increases from lower 
values to higher values. Convergence study of the numerical 
solution can be made by comparison with the known exact solution 
for the limiting. case of straight rectangular channels K=0. By 
using the mesh size of M x N = 20 x 10 for curved square channels, 


the values for (fRe), and (fRe), are found to agree with each 


2 
other to within one per cent up to K=202. For higher Dean 
numbers the mesh size of 20 x 10 is found to be inadequate and 
the mesh size of 40 x 20 is used. The numerical results are 
listed in Tables 1 to 4 for the aspect ratios y = 1,2,5 and 0.5, 
respectively. The difference in values based on two alternative 
methods of computing fRe can be traced back to the difficulty in 
calculating the accurate velocity gradient at wall, (aw/an), 


using the three-point finite difference approximation, because 


of rather steep velocity gradient near the wall. After 
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considerable numerical experiments, it is established that the 
values for (fRe), based on overall force balance are reliable 
and believed to be accurate within two per cent for higher Dean 
number range. It is also found that the curvature effect in the 
definition for (fRe), in equation (20) is rather small and can 
be neglected practically. However, the value for (fRe), in 
equation (20) is found to be closer to (fRe). than that using 
equation (20a). Based on the above discussion, the numerical 
values from (fRe), are taken to be the solution in the present 
analysis. The computing time required is less than two minutes 


on IBM 360/67 for each Dean number investigated. 


2.3 Results and Discussion 
2.3.1 The effect of Dean number on secondary flow patterns 

The effect of Dean number on secondary flow pattern 
in the form of stream function contours is shown in Fig. 2 for 
fully developed flow in a curved square channel. The occurrence 
of two vortices for curved pipes or rectangular channels is well 
understood. However, Fig. 2(b) shows that at K=202, an additional 
counter-rotating pair of vortices appears near the central part 
of the outer wall. This situation was also pointed out in 
earlier work [11]. In interpreting the graphical results, it is 
well to note that the values of stream function at the eye of the 


vortice represents the strength of the secondary flow. The 
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intensity of the secondary vortices shown in Fig. 2(b), (c) and 
(d) is seen to be comparable to that of primary vortices. It is 
of particular interest to note that at K=520, the additional 
secondary vortices disappear completely as shown in Fig. 2(e) 
and the secondary flow pattern also changes appreciably. The 
appearance of the secondary vortices near the central outer wall 
can be understood since the flow field is known to be unstable 
there due to centrifugal instability [13]. Apparently, the 
centrifugal instability problem in curved rectangular channels 
needs further investigation. The shift of the location of 
vortex center with the increase of the Dean number is also of 


interest. 


AS PVCU ie effect of aspect ratio on secondary flow patterns 

The secondary flow patterns for the aspect ratios 
¥ 2.2.5 and 0.5 are shown in) fig: 3, 4. and 9 respectively, for 
various Dean numbers. In presenting the graphical results, 
an attempt was made to illustrate the new aspects of the secondary 
flow patterns and the circumstances under which the counter- 
rotating secondary vortices appear in addition to the primary 
vortices. With aspect ratio y = 5, a pair of secondary vortices 
with weak intensity appear at a rather low Dean number K = 76. 
With y = 5, one clearly sees that the eye of the primary vortices 
moves toward the upper or lower wall with the increase of Dean 


number. 
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In order to gain further insight into secondary flow 
patterns, the distributions of the secondary flow velocity 
components: usv. ‘along x = "50. 35, 04/025 and“y =" 0.1, O25 ina 
Square channel are shown in Figs. 6 and 7 for Dean numbers K = 100 
and 488, respectively. It is useful to note that the solid 
line represents the distribution for u-component and the dashed 
line for v-component. 

The direction of the secondary velocity component is 
self-evident and is not indicated in the figure. Based on the 
results shown in Figs. 6 and 7 and the existence of the additional 
counter-rotating vortex rolls, the validity of boundary-layer 
approximation for high Dean numbers flow regime needs re- 
examination. For example, at K = 488 the reverse flow appears 
for u and v near the central outer wall due to the additional 


counter-rotating vortices. 


2.3.3 Axial velocity distributions and comparison with experimental 
data 
The axial velocity profiles along the horizontal and 
vertical center-line are shown in Figs. 8 to 12 for the aspect 
ratios y = 1,2,5 and 0.5 with Dean number as parameter. Experimental] 
data [10, 14] are available for the case of curved square channel 
and some representative comparisons are shown in Figs. 8 and 9. 


Some discrepance is observed between theoretical and experimental 
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results but it would be rather difficult to explain the exact 
reasons. For the case of K = 368 shown in Fig. 9, an additional 
pair of counter-rotating vortices already appears and the probe 
for velocity measurement may disturb the flow pattern there. No 
attempt was made to compare the present numerical results with 
other experimental velocity distributions shown in Figs. 7 and 8 
of [14] since the Dean number reported there can be obtained 
only by trial and error method. However, a qualitative comparison 
is possible for some Dean numbers. For example, the vertical 
velocity profile for K = 488 shown in Fig. 8(b) agrees well with 
the profiles for K > 251 shown in Fig. 7 of [14]. The behavior of 
horizontal central velocity profiles for K = 202 and 488 is con- 
sistent with that of stream-function contours shown in Fig. 2(b) 
and (d), respectively. 

In Fig. 10, the effect of Dean number on axial 
velocity profiles is seen to be rather regular indicating rather 
Stable flow field up to K = 278 for y = 2. At y = 5, the effect 
of Dean number on axial velocity profile is seen to be rather 
weak up to K = 123 investigated. For y = 0.5, the behavior of 
axial velocity profile is qualitatively similar to that of square 
channel. For K>176, the location of the maximum axial velocity 
moves progressively toward x = 0 and its magnitude decreases as 
Dean number increases. Although not presented here, it is noted 
that the horizontal central velocity profiles for K = 520 and 


618 become quite irregular for y = 1 and it is tempting to consider 
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this phenomenon as the transition from laminar to turbulent flow. 
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The friction factor results are summarized in Fig. 13 
together with the available theoretical and experimental results 
[9, 10, 11 15 ], for square channels only, plotted for comparison. 
In Fig. 13, the aspect ratio effect on friction factor results 
is evident. The difference in the results for y = 0.5 and 2 is 
of some interest. The reason can be explained from the axial 
velocity distributions. For y = 1 and K<200, the present numerical 
result agrees very well with earlier results [11]. The effect 
of curvature ratio on friction in square channels is also shown 
in the inset of Fig. 13 using somewhat magnified scale. The | 
trend of the curvature ratio effect on friction factor result 
is similar to that observed in curved pipes [5, 16]. The 
difference between friction factor results for r = 10 and 100 
is rather small. In view of the curvature ratio effect noted in 
[17] for curved tubes, one may conclude that when r>10, Dean's 
original formulation [1,2], neglecting curvature ratio effect, 


is practically applicable also for curved rectangular channels. 


2.3.5 Correlation equation for friction factor 
It is expected that the friction factor results for 


curved rectangular channels can be correlated as a power series 
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in K as demonstrated for the case of fully developed laminar 

flow in curved pipes [18]. The correlation equation is useful 
for design purposes and the results for aspect ratios y = 1, 2, 5 
and 0.5 based on the present numerical results and using the 


least Square method are given as, 


By a1 fet, SA0p8~1 ~3/2 -2 
Periee a qh (1 + C, K +C, K+ Cy K + Cy K ) 
where (40) 
Y Cy C, Co C, Cy 
] 0.1278 -0.257 0.669 187.7 -512.2 
Z O22736 -24,79 Sco.c -159] 2728 
5 0.0805 ~5.218 104.4 -202.8 0 
0.5 0.0974 4.366 -13.56 Sis -182.6 


The maximum deviation is less than two per cent for all cases in 


the available Dean number ranges as shown in Fig. 13. 


For the case of curved circular pipes, the friction 


factor ratio for the same Reynold's number in curved and straight 


pipes, cite » is known to depend on curvature ratio [5, 16]. 


However, the curvature ratio effect on friction factor results 


is not expected to be detected readily by experimental invest- 


igation [15]. In interpreting the predictions based on correlation 


equation (40), it is well to note that the correlation at higher 


Dean numbers, for example, for y 


results for r = 


1, is baSed on numerical 


10 and 4 listed in Table 1. 
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correlated equation from Ludwieg's experimental data [9] shown in 
Fig. 13 may be attributed to the curvature ratio effect. The 
experimental data shown in Fig. 13 for the range K = 10° % 10° 
are taken from [15] for r = 17.5 only and the agreement there 

is seen to be reasonable. In this study, no attempt was made 

to develop correlation equation considering curvature ratio 
effect. In high Dean aumberrange, the slope of the curve from 
the present correlation seems to agree well with that of the exp- 


1/2 


erimental correlation equation fel ihe = 0.107K proposed by 


Baylis [15] for square channel. 


~ 2.4 Concluding Remarks 

L. Fully developed laminar flow in curved rectangular 
channels is approached by numerical method considering curvature 
ratio effect in the formulation. It is found that when r>l0, 

the curvature ratio effect can be neglected practically and Dean's 
original formulation is applicable. 

ran The present analysis based on the complete Navier-Stokes 
equations reveals that an additional pair of secondary vortices 
appears near the central outer region of the channel at a certain 
higher Dean number depending on the aspect ratio. This phenomenon 
is consistent with Dean's instability problem for curved parallel- 
plate channel flow. The existence of two pairs of primary and 


secondary vortices poses some questions regarding the adequacy of 
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the boundary layer approximaticn for the high Dean number flow 
regime. Within the scope of present investigation, the secondary 
vortices disappear at K = 520, for example, for the curved 

Square channel. 

Se The numerical results are tabulated for future reference 
and the correlation equations for the friction factor are developed. 
The graphical results are designed to illustrate new aspects of 
flow fields in curved rectangular channels. For square channels, 
the numerical solution covers the Dean number range 5 ~ 715. 

It is noted that a convergent solution can be obtained by the 
present numerical method for higher Dean numbers using finer mesh 
sizes but at the expense of considerable computing time. 

4, The solution of the to problem in curved rectangular 
channels is basic to various heat and mass transfer problems. 

Some related thermal entrance region problems will be studied 


in the next three chapters. 
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CHAPTER ITI 


GRAETZ PROBLEM IN CURVED SQUARE CHANNELS 
3.1 Formulation of the Problem 


The classical Graetz problem is the theoretical analysis 
of laminar heat transfer in the thermal entrance region of 
Straight tubes or channels. This analysis has been extended 
by several authors to different wall conditions and various real 
effects which were not considered originally by Graetz. The 
Graetz problem in curved square channels is characterized by the 
inclusion of the secondary flow effect, induced by the centrifugal 


force, in addition to the axial main flow effect. 


The following assumptions are made in order to limit the 
scope of the present analysis: 

1. The physical properties are constant and free 
convection effect is negligible. 

2. Viscous dissipation effect is negligible. 

3. Heat sources do not exist. 

4, The axial conduction term is negligible for large 


Peclet number, (Pe>50). 


With assumption 1, the energy equation is seen to be 
uncoupled with the momentum equations. The flow problem solved 


in Chapter II can be applied to the present problem as well. One 
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can just concentrate on the discussion of energy equation. 

The governing equation for the energy transport in curved square chan- 
nels is given below using the same coordinate system as that used 

in Chapter II for flow problem. The derivation of the energy 


equation is given in Appendix A. 
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(41) 


The thermal boundary conditions to be considered in this study 
are: 


(1) Constant wall temperature 


T= M2 at entrance (z<0) 

Pels at wall 

af = 0 along the line of symmetry y = 0 
(2) Uniform wall heat flux 

Laka at entrance (z<0) (42) 

Pe -k(St), = const. at walls 

af = 0 along the line of symmetry y = 0 
(3) Uniform convective boundary condition 

Lila at entrance (z<0) 

-k(2),, = TA Bret at walls 

aT _ g li f t = 0 

eile along the line of symmetry y = 
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Introducing the following non-dimensional quantities: 


X = De x, Y = De y, Z = De Pr Re z, R= Der 
= oe ie Ur 
U = De Us V De Ve W = Ww 
1-139 (Tea )e for constant wall temperature case, 


al. = (q, De/k) 6 for uniform wall heat flux, 


and oie = (T)-T,)0 for wall convection, 


where 
De = 2abfatb), Re = W De/v 


equation (41) becomes 


2 2 
00 30 W Pee Oe at \) ps eee ems 1 06 
Uo oy It+x/r Pr 3z Pr sea eore rUlax/ Py) 9x J 
xX dy 
(43) 
and the corresponding boundary condition becomes 
(1) Constant wall temperature case 
6 = ] at entrance (z<0) 
6 = 0 at walls 
ae 0 along the line of symmetry y=0 
(11) Uniform wall heat flux case 
06=0 at entrance (z<0) 
oo8-s . 
pet ] at walls (44) 
o8 0 along the line of symmetry y=0 
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(III) Uniform convective boundary condition 


e= ] at entrance (z<0) 
: cS Ags t walls, where Bi=h De/k 
7 ‘i at walls, where Bi e 
foe 0 along the line of symmetry y=0 


The Nusselt number variation in the thermal entrance 
region is of considerable importance in the present study. The 
Nusselt number Nu = h De/k can be obtained in two ways for each 
boundary condition as: 

(a) Constant wall temperature case 


Considering the wall temperature gradient, one obtains 


(25) (45) 


Using the dimensionless quantities introduced above, one obtains 
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Following the usual definition for the Nusselt number, one has 


ee 

on 
hDe _ W 

Nu = =O 


48 
] 64, ( ) 


where 6, is the bulk temperature which is defined as: 


Afi woot 
ike 


On the other hand, an energy balance of bulk flow with axial 


length dZ yields, 


o Cp THUR a= at A dZ=h (T,-T,) dZs 
W oT 
fe p Cpl x 
h=¢ aS ; ow 
W b 
where 
W of os 
( THR az ? “K Sf THX/R a7 
Using the same dimensionless quantities, one obtains 
te Wen = 36> 
thes P RePrde 1tx/T Oz 
Lees La oe, 
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Hence, one has 
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(b) Uniform wall heat flux case 


Performing the same procedures as for case (a), one obtains 


Nu, = = (53) 
~6 
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cs ee _ WwW 88 : 
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(c) Uniform convective boundary condition 


Following a similar procedure as noted in case (a), one obtains 
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Simpson's integration formula is employed to carry 
out the calculations of the average quantities appearing in 
equations (48) to (56). Simple backward difference approximation 
is used to calculate the temperature gradient, 96/az. As to 
the wall temperature gradient calculation in equation (48), 


three-point finite-difference approximation is used. 


The bulk temperature development relationship in 
the uniform wall heat flux case is of special interest. 
Considering the energy balance for the axial length dZ, of the 


channel, one has 


o] 


a ) 
ep © TGR a7 A dZ = q.+S+dd 


cas oe 
THOR 92 pcp A 


After using dimensionless variables, the above equation reduces to 


= 4 (58) 
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Equation (59) can be written in the following way: 


or 


From the non-dimensional quantities, one has the following 


relationship: 


Luann =] (61) 


Equation (60) is divided by the left hand side of equation (61), 


then one obtains 


06 
eae 
Favret 


Two separate methods of calculating Nu and equation (62) provide 


a means of assessing the accuracy of the numerical solution. 


3.2 Finite Difference Solution 


The velocity components, u,v, and w are available from 
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the numerical solution of the fully developed flow problem 

in Chapter II. Consequently, the finite-difference grid-size 

MxN used for the energy equation is the same as that used in flow 
problem. With u=v=0 and 1/r=0, the present Graetz problem 

is seen to be analogous to the unsteady state two dimensional 
heat conduction problem. However, the convective terms due 

to secondary flow becomes increasingly important in equation (43) 
as Dean or Prandtl] number increases. The alternating direction 
implicit (AD1) method due to Peaceman and Rachford [22] is 
employed to transform equation (43) into a set of algebraic 


equations. 


The numerical procedure consists of two steps for a 
given axial step distance Az. The entire temperature field is 
known at the previous axial step z (nth step). The half axial 
step computes the temperature values for all grid points at axial 
distance z + 0.5 Az, just considering the variation in the x 
direction. The full axial step computes the final temperature 
at the fully advanced axial distance z + Az, proceeding another 
0.5 Az after half axial step, by using the half axial step values 
and considering only the variation in the y direction. Using 
second-order central difference approximation for all derivatives, 


one has the following finite difference linear algebraic equations: 
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(1) Half axial step (implicit in x) 

Indicating the known temperature at the nth axial 
distance (z) by 6” and the unknown temperature at half axial 
distance (z + 0.5 Az) by a single asterisk, the finite difference 


analog of equation (43) can be written as: 
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The above equation can be rearranged in the following form: 
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ue can be obtained by solving the above equation using 
Gaussian elimination method [12]. With ee available, then one 
proceeds to the next step. 

(2) Full axial step (implicit in y) 


* kk 
With F j as the known values and F j as still unknown 
b $ 


values, the finite difference analog equation becomes 
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Again, using Gaussian elimination method, one finishes a complete 
axial step calculation. Therefore, the (n+1)th step temperature 
distribution is obtained. Repeating these two steps, the 


temperature values in the thermal entrance region can be obtained. 


The corresponding boundary conditions in finite difference 
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form are: 


(a) Constant wall temperature case 


Ore TOM) ack hs 
= 0 for i=1,...M+] (67) 
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(b) Uniform wall heat flux case 
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(c) Uniform convective boundary condition 
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It must be pointed out that in the calculation of 
wall temperature for both uniform wall heat flux and wall 
convection cases, the governing equation takes the following 
simple form: 


36 3¢ 


= ] " 
peciatye cls a) Texte xe (70) 


In finite difference analog, it becomes 


1 ] 2 2 
( i )e.,, «-( —z t+ —s )e. . + 
Axe 2r (14x. /1) Ax ee Axe ay? Tigh 
( ee 1 ie 
Axe 2r(1+x./r) i-1,j 
+o (05 jay? 45-1) = (71) 
Ay sJ oJ 


By applying the boundary conditions at wall, the wall temperature 
can be computed after the temperatures in the interior region are 


found. 


3.3 Convergence of the Numerical Solutions 
The stability of the numerical solution depends on grid- 
size (MxN) and the axial step size (Az). For the flow problem in 


curved square channel, the grid-size of 20x10 and occasionally 


in | m4 5 
wa r 
Ray 
ae a 5) 7 


7 : CA : + aaa Ru he . 
Yo notis!yolso af oe pau UY TEN 


a Cs 


i iy pees ; 


— up Lh! a 

r o7t tr: Ae 

oe Fad ‘ x ia Bp iy vo 
ffew bas xuf ‘$99 | ba 5 baie. 


be . 


= mH ok ean . 
entwilfo? ald eotat a os ork 


van 


- a. 
9% oa #4 34 ; 


~~ 


o 
ry 


Pg 
Fe" 
ay 


> : xe : ie 1\x# = - - : 
- ae : : 
} i ° 


Reet “aie 


zomozsd $f .gotans 8008 na te 


oes 


ee ara” 


oy 1} 3 lew 38 4+ TDTOS Shae 4 
; f ; pit fei Th 
ft 
motiuie2 f42 isi 
shraqeb aoryufee | extvomun SAS to 
rt - ae a, ss > are 
ati » wilt oat YOU .(Sa) este qeze 


vilsnotesose bas OFXOS to estre-b ri ‘ond? . fons 28 a1 6up 2 


abt 


44 


40x20 was found to be satisfactory from the viewpoint of accuracy 
and computing time. For the present problem, the following scheme 
of axial step variation obtained by trial and error method leads 


to stable and convergent solution. 


Axial step, n Step size, Az 
1 ~ 100 5x10°> 

101 ~ 450 1074 

451 ~ 650 2x10" 

651 ~ 950 5x10" 
n> 851 iow 


The solution is terminated when the asymptotic 


solution which will be discussed later, is reached. 


The accuracy of the numerical result can be assessed 
by noting the known asymptotic Nusslet number of 2.98 and 3.09] 
for the constant wall temperature and uniform wall heat flux 
cases, respectively. The present calculations agree very well 
with the foregoing values. For the case of constant wal] 
temperature, the difference in the values for Nu, and Nuy from 
equations (48) and (52) is found to be relatively large when 
Prandtl] or Dean number is large. For example, the asymptotic 
Nusselt number Nu, and Nu. deviates from the average value 
by 7 per cent with Pr=0.7, K=368 and MxN=40x20. Similar 


observations is also noted in [16] for the case of curved pipes. 
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By further reduction in mesh size and axial step, it is confirmed 
that the values from Nu, iS a reasonably accurate one and 
consequently the numerical value from Nu, is used in presenting 
the result. It is suspected that the numerical evaluation of 
06/az 1S a source of difficulty for Nu, in the constant wal] 
temperature case. For the case of uniform wall heat flux, Nu, 
and Nu, are seen to be identical if the relationship represented 
by equation (62) is used. It is also clear that the numerical 
value from Nu, 1S more accurate than Nu,. In contrast to the 
case of constant wall temperature, the two values from Nu, and 
Nu, are found to agree with each other very well. At a higher 
Dean number K=488, the deviation from the average value is 

found to be 3 per cent with Pr=0.1 and MxN=20x10. The numerical 
result also checks well with equation (62). For the case of wall 
convection, Nu, and Nu. agree with each other just as good as in 


the case of uniform wall heat flux for all parameters investigated. 


From the foregoing discussion, it is concluded that 
the numerical values from Nu, are sufficiently accurate for these 
thermal boundary conditions and the accuracy is considered to be 
within about 2 per cent. The computing time required in obtaining 
a complete result for each Dean number is approximately 3 minutes 


with MxN=20x10 and 9 minutes with MxN=40x20 on IBM 360/67 system. 
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CHAPTER IV 


RESULTS FOR CONSTANT WALL TEMPERATURE 
AND UNIFORM WALL HEAT FLUX BOUNDARY CONDITIONS 


4.1 Results for Constant Wall Temperature Case 


4.1.1 Temperature field development 


In this study, the effect of Dean number on the 
development of temperature field in the thermal entrance region 
is of special interest. The temperature profile developments 
along the central vertical and horizontal axes are shown in Fig. 14(a) 
to (d) for Dean numbers K = 14, 100, 368, and 488, respectively, 
with Pr = 0.7. In examining the graphical results, one should 
note the gradual distortion of the profile and the shift of the 
location of maximum value for 6 as the Dean number increases. It 
is pointed out that at K = 368 and 488 two pairs of counter-rotating 
vortex rolls already appear and the detailed secondary flow 
patterns are shown in Chapter II. The development of isothermals 
at several axial locations for Pr = 5 and K = 100 is illustrated 
in Fig. 15. The development of the plume-like behavior in the 
form of warm current penetrating into the cold fluid core region 
and the appearance of the two eyes for isothermals are note- 
worthy. The phenomenon is similar to that observed in curved 


circular pipes [19, 20, 21]. The effect of Dean number on bulk 
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temperature distribution along downstream direction is shown in 


hig. 6ator Py = 0.7 sand 10.05. 


4.1.2 Heat transfer result 


The variation of local Nusselt number in the thermal 
entrance region is of practical importance and the results for 
Peaa?07 hes 02/5 Seand 500;are%shown- inFigs. 17 andl8.. AtePr = 0.1; 
the general feature of local Nusselt number variation is similar 
to that of the classical Graetz solution (K=0). The behavior for 
Nu can be explained from equation (43). At small Prandtl] number 
the axial convective term dominates over the convective terms 
due to secondary flow. Consequently, the curves for different Dean 
numbers are nearly parallel and straight in the Leveque solution 
region near the thermal entrance (z=0) and the difference in 
Nusselt number values is caused by the different degree of axial 
velocity profile distortion for w due to Dean number effect. For 
small Prandtl numbers, the role of axial convective term in the 


energy equation is important. 


At Py ="077, one can already see the existence of 
minimum Nusselt number values at some downstream location depend- 
ing on Dean number. At this particular location, apparently, 
the axial convective term and the convective terms due to 
secondary flow are equally important and one may conclude that 


the entrance effect and the secondary flow effect are in balance. 
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After passing the minimum Nusselt number point, the secondary 

flow effect dominates over the entrance effect and the local 

Nusselt number continues to increase until the local maximum 

values is reached. Eventually, the Nusselt number values remain 
constant further downstream and one may conclude that the limiting 
asympotic value, corresponding to the fully developed temperature field 
is reached. It is clear that the thermal entrance length 

decreases with the increase of Dean number. For the case of 

Pr = 0.7, it iS instructive to contrast the Nusselt number 

variation with temperature profile development shown in Fig. 14 


and the bulk. temperature distribution shown in Fig. 16. 


A characteristically different kind of Nusselt number 
behavior appears when Prandtl number is large. For example, at 
Pr = 500 and a relatively low Dean number K = 14, the secondary 
flow is rather weak and the Nusselt number variation in the 
Leveque solution region is identical with the case with K = 0. 
However, at some downstream position the secondary effect appears 
and the deviation occurs. After reaching the local maximum value 
for Nu, a fluctuating phenomenon appears before a stationary 
value is reached further downstream. The present fluctuating 
Nusselt number phenomenon is similar to that reported in the 
earlier numerical studies [16, 19] for the case of curved pipes. 
It should be pointed out that numerical difficulty was encountered 
for K>14 with Pr = 500. From the foregoing discussion, one can 


also understand the Nusselt number behavior for Pr = 5 as shown 
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invhid, 18. 


4.2 Results for Uniform Wall Heat Flux Case 


Temperature field developments and heat transfer results 
Similar to those presented earlier are shown in Fig. 19 to 23. 
The behavior of temperature profile development shown in Fig. 19(b) 
can be understood better if it is contrasted to the corresponding 
isothermal pattern shown in Fig. 20 and the distributions for oy 
and 6, shown in Fig. 21. For the case of uniform wall heat flux, 
apparently the left-hand corner is a hot spot and the wall 
temperature is highest there. When the temperature field 
becomes fully developed, the axial distributions for 8h and 6 
become parallel to each other and correspondingly the local 
Nusselt number approaches a constant value depending on Prandtl 
and Dean numbers. The thermal entrance length is of practical 
interest and can be determined from Figs. 21 to 23 by noting the 
above fact. The local Nusselt number result shows that for a 
given Prandt] number, the thermal entrance length decrease with 
the increase of Dean number and the Prandtl number effect is also 
seen to be similar for a given Dean number. The roles of Dean 
and Prandtl numbers on heat transfer result are seen to be similar 
and this observation can also be seen from equation (43). The 
relative magnitude of the convective terms due to secondary flow, 

90 30 


Pr(u erm 1 oy ), increase with the increase of Prandtl] or Dean 


number. 
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The behavior of local Nusselt number in the Leveque 
solution region and the occurrence of minimum and subsequent local 
maximum Nusselt number can be explained. However, the fluctuating 
Nusselt number phenomenon before approaching the asymptotic value 
as shown in Fig. 23 cannot be explained readily. The reported 
numerical results [16, 19] and the present work show qualitative 
agreement regarding the fluctuating Nusseit number phenomenon. 
Apparently, the fluctuating phenomenon is caused by the convective 
terms due to secondary flow in the energy equation (43). Consider- 
ing the fact that numerical difficulty also arises when Prandtl 
or Dean number is very large, one cannot rule out entirely the 
possibility that the phenomenon is caused by numerical instability 
[20, 21]. At present it is not clear why the Nusselt number 
decrease faster along the axial direction than the case of 
Straight channel. Perhaps the question can only be answered by 
very precise and accurate experimental measurements. However, 
the fluctuating phenomenon is probably of theoretical interest 
only and in practical applications the asymptotic Nusselt number 


is important when Prandtl or Dean number is large. 


- 4,3 Asymptotic Nusselt Number Results 


The asymptotic Nusselt number is of practical import- 
ance and the numerical results for the two thermal boundary 


conditions are listed in Table 5 and plotted in Fig. 24. The 
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experimental data [14] for the case of uniform wall heat flux 

with Pr = 0.71 are also shown there for comparison. The agreement 
with the present analysis is seen to be reasonable. This confirms 
the good accuracy of the numerical analysis noted earlier. The 
difference between the two thermal boundary conditions f seen to 
be relatively small for Pr = 0.7 and 5. In view of this fact, 

the numerical results for the two boundary conditions were 
correlated by the following equation within a maximum deviation 


of 12% using the least-square method. 


Nu = 0.152 + 0.627(K!/* prl/4) | o.7<Pr<c5, 14<K<500 


(72) 
V2) 


The parameter (KPr is based on earlier study [23] for curved 


pipes. 


4,4 Concluding Remarks 


1. The present analysis includes the curvature ratio effect in 

the formulation. However, recent work [24] for fully developed 
laminar forced convection in curved tubes with the thermal boundary 
condition of axially uniform wall heat flux and peripherally 
uniform wall temperature concludes that the curvature ratio 
parameter has a negligible effect on the average Nusslet number 

in the range of R/a = 10 to 100. For the present problem involving 


curved square channels, the curvature ratio effect is also expected 
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to be negligible for Nu with r>10. 

2. In examining the graphical results presented, it is well to note 
that two pairs of vortex rolls already appear at Dean number K = 202, 
for example. In the results of Chapter II, no attempt was made 

to determine exactly the Dean number corresponding to the onset of 
centrifugal instability. 

3. The present numerical results reveal fluctuating Nusselt number 
behavior before reaching the asymptotic value when the Prandtl 

or Dean number is large. The phenomenon is similar to that observed 
in curved circular pipes [16,19]. An attempt was made to clarify 
the heat transfer mechanism particularly the Prandtl and Dean 

number effects in the thermal entrance region of curved square 
channels. 

4. In order to limit the scope of the present investigation, the 
aspect ratio effect on the present convective heat transfer problem 
in curved rectangular channels was not studied here. However, the 
analysis can be made by using the flow fields reported in Chapter II. 
5. Some local Nusselt number results in the thermal entrance 

region are tabulated in Appendix C.1 and C.2 for possible future 


reference. 
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CHAPTER V 


RESULTS FOR UNIFORM WALL CONVECTIVE 
BOUNDARY CONDITION 


5.1 Temperature Field Development 


The development of the temperature field along the 
central and vertical planes of a curved square channel is shown 
in Fig. 26 for Pr = 10.05, K = 368 and Bi = 1, 10. The 
corresponding wall temperature distribution oe around the channel 
periphery is shown in Fig. 27. It is noted that at K = 368, 
an additional pair of vortex rolls already appears near the 
central outer region of the channel. It will be instructive 
to examine the temperature profile development with local 
Nusselt number variation which will be shown later. Because 
of the nature of present convective boundary condition, the 
developing wall temperature distribution for 6, is of special 
interest. It is clearly seen that the value of e, is lowest at 
the inside corners. The temperature values are lower near the 
corners and the corners points may be regarded as cold spot. 

It should be pointed out that the secondary flow is rather weak 
near the corners. Noting the relationship 96/on = - Bi e., from 
equation (44) and the fact that the Biot number is constant, it 


is seen that the distribution of the wall temperature gradient 
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a6/on is also similar to the distribution for Os 


The axial variation of bulk temperature 8, for, Pr = 
0.03, 0.7 and 10.05 are shown in Figs. 28 to 30 with Dean 
number as parameter. At Bi = 1, the Dean number effect is 
relatively small. These figures also illustrate the Biot 
number effect. With perfect insulation (Bi>0), the bulk 
temperature 8b remains to be one. For a given Dean number, 


the Prandti number effect is seen to be rather small. 


5.2 Local Nusselt Number Result 


The local Nusselt number variation in the thermal 
entrance region of curved square channels are illustrated in 
Figs. 31 and 32 with Dean number as parameter for Pr = 0.03, 0.7, 
and 10, respectively. As shown in Fig. 31, the Biot number 
effect on local Nusselt number is rather small and at K = 368, 
for example, the values with Bi = 1 is only slightly higher than 
those for Bi = 10. Although not shown in Fig. 31, the value 
for the limiting case (Bi>~) is only slightly lower than the 


value for Bi = 10. 


At low Prandtl number such as Pr = 0.03, the local 
Nusselt number behavior is similar to that of straight channel 
K = 0. This phenomenon has been discussed in Chapter IV for 
constant wall temperature and uniform wall heat flux boundary 


conditions. For Prandtl number 0.7 and 10.05, the local Nusselt 
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number fluctuates before the fully developed stationary value is 
reached. This fluctuating behavior has already been pointed out 
in Chapter IV also. The asymptotic Nusselt number is of practical 
interest and the results are given in Table 6. It is noted 

that the Biot number effect is rather small and a maximum 
deviation from the uniform wall temperature case (Bi-~) is 


only 4 per cent for a given Dean number. 


5.3 The Effect of Secondary Flow on Liquid Solidification-Free 


Zone 


In the case of uniform convective cooling at the outer 
surface with ambient temperature T. lower than the liquid freezing 
temperature Te and in view of the wall temperature distribution 
snown in Frg. (275 bt. 1s cleay that liquid solidification will 
start at the inner corner point B when the freezing temperature 
is reached there at some downstream distance referred to as the 
length of solidification-free zone or the solid-phase entrance 
(Zz = Ze). When the solid-phase surrounds the channel inner wall 
completely, the thermal boundary condition becomes isothermal. 
The solidification of liquid metal or water in curved channels 
is of considerable practical importance. However, in this 
investigation only the prediction of solidification-free zone 


will be considered. 
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For a study of solidification-free zone, it is 
convenient to introduce the superheat ratio parameter [28] 


defined as, © = (T-Te)/(Te-T,,) = - 1 + 67/, where Opel. It 


f 
is noted that solidification begins when the wall temperature 
§, at inner corners (cold spot) equals 6-. Consequently, the 
Superheat ratio « corresponding to the solidification-free zone 
length Z~ can be computed readily. The results are plotted in 
Fig. 33 for Pr = 0.03 (liquid metal) and 10.05 (water) with 
Dean number as parameter for the cases Bi = 1 and 10. Since 
the present analysis is based on the constant physical property 
assumptions, with Pr = 10.05 for water, for example, the 
temperature difference (T,-Te) is practically limited to about 
20°C. Furthermore, with T, = -0.5°C and -100°C, « becomes 40 
and 0.2, respectively. In interpreting the results shown in 
Fig. 33, it is well to note that as e©« approaches zero and 
infinity, Ze also approaches zero and infinity, respectively. 
The trend for K = 0 is similar to that reported in [28, 30] for 
Straight tubes. With Pr = 0.03 and Bi = 1, the secondary flow 
effect seems to be small. It is sufficient here to use the 
case of Pr = 10.05 and Bi = 10 to explain the secondary flow 
effect on solidification-free zone. Apparently, the character 
of Dean number effect changes at the threshold value of e=37. 
Noting that the value of e« decreases as the ambient temperature 


ue decreases, with a given e < 37, the secondary flow tends to 


‘at at anos. 99 te tte ryfott: 
[8S] ves seie199 atts ant sits ay sou 


ay : 


t ,f>,8 svariw ee kes 
i ( si 7 
vutsvoqnat (iow ett oadantpad 0h 143 tte 
Payee . ne 
aid yfingupshnod 40 ee ‘8 
; iti ioe if o 
' ttol Sfyeoy oT 
footew) 20.0% bas (Te :btuptt) 20,0898 bad 
20250 ort “107 sat ome: ia ies 1 Sedinae 
7 bey ut Pie 
vra 27ezte2ne2 Sf ya bsend 2t atey s 209 ¥ 
9 20,0n 4 4 Aiw ie ok 
a 
a4 nal ‘ef vifsor3 ard ir 31% 1) pat 8 r¥ 
) 1 | : ; = 


i ANS ++ 
) : 
~ D> P “ 
7! [ c J 
4 
ious? 2 3 3 » find 
ee r = tf } i 
~ | 3 
At 
‘ fl ra +, 4 { 
+ 
2 7 iA ate ST 
: Vii } : < wiih 


tnatdns ah? 28 292s 79eb 3 Fo i tae | dent 


+ \vrsbroase ont NE > o navig & ast a 2 


= 
yy 


o/ 


mix the warmer liquid and correspondingly, the solidification- 
free length increases with the increase of Dean number. On the 
other hand, with e > 37 (or relatively higher TH) the secondary 
flow tends to bring the cold spot temperature to the freezing 
temperature at shorter entrance distance Zee Further physical 
insight can be gained by considering a given value of Z_ and 


the corresponding different values for different Dean numbers. 


5.4 Aspect Ratio Effect 


The results for the typical local Nusselt number and 
solidification-free zone for Pr = 10.05 (water) and 0.03 (liquid 
metal) are also presented in Figs. 34 and 35, respectively, for 
aspeciiratiosuye= 11 .e2shorahdhOv5qwi thiBie=tieandedOsealthe 
Nusselt number behavior is characteristically different between 
Pr = 0.03 and 10.05. The fluctuating Nusselt number phenomenon 
before approaching the asymptotic value for Pr = 10.05 and 
K>171 is somewhat similar to the earlier investigation [16, 19, 31] 
but the number of fluctuating cycles before the complete damping 
is far fewer than that shown in Fig. 9 of [16] for the constant 
wall temperature case with Pr = 10 in helical coils. The aspect 


ratio effect can only be seen clearly with comparable Dean number. 


5.5 Concluding Remarks 


1. The Graetz problem in curved rectangular channels 
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without the assumption of large curvature ratio R/De is solved 
numerically for the uniform convective boundary condition. The 
solution also yields the length of liquid solidification-free 
zone. The numerical result reveals that solidification will start 
at the inner corners which represent a cold spot. One may also 
reason that solidification will start at inner wall for the 


practically important case of curved circular tubes. 


2. The Nusselt number behavior in the thermal entrance 
region is characteristically different between large and smal] 
Prandti numbers. For given values of Prandtl and Dean numbers, 
the Biot number effect on local Nusselt number is found to be 
rather small within the range of present investigation. The 
fluctuating Nusselt number behavior before the eventual damping 
of the cyclic behavior and the approach of the asymptotic value 
is similar to that reported in earlier investigations. Both the 
Prandtl and Dean numbers have the effect of decreasing the 


thermal entrance length. 


3. For a given value of Biot number, the secondary flow 
effect on the length of liquid solidification-free zone is greater 
for larger Prandtl] number. For given values of Prandtl! number and 
superheat ratio, the Dean number effect on solidification-free 
length is completely opposite depending on whether the solidification 
occurs in thermally developing or fully developed zones. The 


present investigation on solidification-free zones also provides 
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some insight into the solidification problem in curved channels 


for liquid metal and water. 


4. Typical local Nusselt number results in the thermal 
entrance region are tabulated in Appendix C.3 for possible future 


reference. 
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CHAPTER VI 


CONCLUSION 


It is found that when r > 10, the curvature ratio 
effect on laminar flow results can be neglected practically. 
In laminar forced convection heat transfer analysis, the 
curvature ratio effect is also included in the formulation. 
However, it is also concluded that the curvature ratio 
effect is expected to be negligible for Nusselt number 


Wrtn Tr -10". 


The appearance of the additional secondary vortex 
rolls at a certain higher Dean number for fully developed 
laminar flow in curved rectangular channels is a new 
physical phenomenon revealed by the present analysis. The 
intensity of the secondary vortices is seen to be comparable 
to that of the primary vortices. Although some experimental 
works on curved channel flows have been reported in the 
literature, the effort appears to have been directed toward 
obtaining experimental data for friction factor and heat 
transfer coefficient only. It is suggested that a flow 
visualization study be conducted to confirm the existence 
of the additional streamwise vortex rolls near the central 


outer wall of the curved rectangular channels. 
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The Prandtl and Dean number effects on the heat transfer 
results are found to be similar for the three different thermal 
boundary conditions investigated. When Prandtl number is 
small, the local Nusselt number variation in the thermal 
entrance region behaves in the same way as the Graetz 
problem in straight channels. When Prandtl] or Dean number 
is large, some cyclic behavior of Nusselt number appears 
before the fully developed asymptotic value is reached. 

For large Prandtl or Dean number case, the fully developed 
region is of practical interest since the thermal entrance 
length decreases with the increase of Prandtl] or Dean 
number. Within the scope of this investigation, the effects 
of different thermal boundary conditions on the asymptotic 
Nusselt number values appear to be small for given Prandtl 
number and Dean numbers. In the case of uniform convective 
boundary condition, the cases of Biot numbers 1 and 10 are 
studied. The difference between the Nusselt number results 
for Bi=] and 10 is found to be small. The limiting cases 
(Bi-0) and (Bie~) correspond to the complete insulation and 


the constant wall temperature cases, respectively. 


It is found that the secondary flow effect on the 
liquid solidification-free zone in the curved rectangular 
channels is significant. For given values of Prandtl 


number and superheat ratio, the secondary flow effect on 
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solidification-free zone is completely opposite depending 
On whether the solidification occurs in the thermally 
developing or fully developed region. Although the 
problem has not been studied for the case of practically 
important case of curved circular tubes, one may expect 
that the trend of the phenomenon would be the same. In 
view of the significance of the secondary flow effect on 
liquid solidification-free zone in curved channels, a 
study of the secondary flow effect due to buoyance forces 


in straight tubes or channels seems to be important. 


Graetz problem in curved tubes has TICKET by 
several investigators in recent years. It is worthwhile 
to point out that hydrodynamic entrance region problem 
in helically coiled pipes has been solved successfully by 
numerical method [31]. This fact immediately suggests that 
the hydrodynamic entrance region problem in curved 
rectangular channels can be approached by numerical method. 
Furthermore, a possibility exists that the entrance region 
problem for simultaneous development of velocity and 


temperature fields can also be solved by numerical method. 
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APPENDIX A 


DERIVATION OF THE BASIC EQUATIONS 
FOR FLOW AND HEAT TRANSFER 
IN A CURVED RECTANGULAR CHANNEL 


The governing conservation equations for a steady, 


incompressible, laminar Newtonian flow in vector form are: 


0 (V . w) = - YP + (VV) (A.2) 


Neglecting viscous dissipation, the energy equation becomes 


0 Cp(V:vT) = kveT (A.3) 


An orthogonal curvilinear coordinate (toroidal) system X,Y,Z, 
as shown in Fig. A.1, is introduced. The relationships between 
the toroidal coordinates (X,Y,Z) and the Cartesian coordinates 


(MaVvez ares 


X = (Ro tek eeose 


7 
ete ie Fee 
ved 
fae 


- ily, | 
| ; a | ar 
< -) 

iv 

ev near 7 vorray 

a 

Ji2WAaT TA a3 ci Os 1 + 

JAMMAHD SASUOMATD IA oi - AW 

- 


t 
' : wis 182 U9 
= P i sQiek. i} 
P VY HY * ' 4 hs I aK 
bal - 
~- 
7, ‘¥ 


¥ A _ 


eee 
) nottsups yp rgsqteeth ati f Rats: 287 O: 


; ie a 
part tv nv) Tsnogortyre 
5 +} ue Pisial i af (EA, ari ow 0) ne 
18) ort brs (Sy ¥.X) 2 SI HATHTOOD febtor MO: 


i, 


114 


yaad 
Ze CR ee etna (A.4) 
¢ = Z/R 


The metric coefficients [32] are: 


1 Ox aX aX 
AAs OR. ay \2 zm 2K_ 
hs ( ay ) + ( aY ) 2 ( JV ) a ] (A S) 


hy = ( SR )h 4 (4 (EY? = (r4xyR)? 


and ce ©. Gs are the unit base vectors in the X,Y, and Z 
directions, respectively. The derivatives of unit base vectors 


are; 


an eas ee (A.6) 
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a ts 
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(b.A) 


Yi r j= , : 
+ He y= 


he ia ie ath 
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(a,A) { = 
5 io ix +f ) a= 


i bits , 1,8 adit nt evotoey cesd thaw eves a) oe a 
eros3ev sesd gFou da 2eviteviveb edT i 


La ic, 


The del operator is defined as: 
de (A.7) 


The Laplacian operation ve = V-V is obtained by the indicated 


repeated application of the operator V. Thus, one obtains 


2 2 
ad a e) ] 1G4lee 


52 
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3 oY hy oZ 3 


The continuity equation is written as: 
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which combined with equation (A.6) becomes: 
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The steady state Eulerian derivative V-V in equation (A.2) is a 


scalar point function and is given as: 
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The quantity V-v V is therefore given as: 
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The pressure gradient VP is: 
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The Laplacian vector vv becomes: 
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From equations (A.13), (A.14) and (A.15), equation (A.2) can be 


written in the following component forms: 
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The equations (A.10), (A.16), (A.17) and (A.18) represent the 
complete set of governing equations for three-dimensional steady 
incompressible flow in curved rectangular channels. For the 
case of fully developed flow, the set of equations can be 
simplified further by noting that the Z-derivatives are zero. 


The governing equations then become: 


Continuity equation 
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X - momentum equation 
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Z - momentum equation 
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The energy equation will be considered next. Using the expression 
V-v from equation (A.11) and the Laplacian operator v from 


equation (A.8), the energy equation (A.3) becomes 
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It is known that the effect of axial conduction term in the energy 
equation can be neglected practically when the value of Peclet 
number is large (Pe>50). Neglecting the axial conduction term, 


the energy equation takes the following form: 
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The Graetz problem (thermal entrance heat transfer) 
in curved rectangular channels is concerned with the solution of 


the above parabolic equation with appropriate boundary conditions. 
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APPENDIX 8B 


COMPUTER PROGRAMS 


He oe oe sk 3 SK SK eK 3 KK ek aK aK KK ok eK KKK KKK KOR KKK KK OK KOK KK KK KK KK OK KK 
LAMINAR FLGW IN CURVED RECTANGULAR CHANNELS 
HK eo oe 9K AK ok ok a a SK a Ko 2 A a KO a I a a aK KC KK eK aK a oe ok ok oe ok Kk 
MsN NUMBER OF GRID POINTS 
R ASPECT RATIO 
0D PRESSURE DROP 
RC RADIUS CURVATURE 
MAXIMUM NUMBER OF ITERATION 
ITMs ISMs IVM MAXIMUM NUMBER OF ITERATIONS 
FR1 PRODUCT OF FRICTION FACTOR AND REYNOLD NUMBER 
1 FGR STRAIGHT CHANNEL 
DIMENSION X(41)sRA1(41)sRAZ2(41) 5S(41541).W( 41 541); 
1VOR(41 541) 
*oS1(41 4.21) »W1(41521)sVOR1 (41.421),U(41.,21) 5V(415,21)>s 
1A1(41)+A2(41) 
*9A3(41)+A4(41)50W(41541) »FI(41) 
DIMENSION S2(41 641) +VX( 41541) 
COMMON /C4/ MN 
COMMON /CS/ MieNl 
REAL MAXS»sMAXT 
READ(53:10) MeNsR 
FORMAT (2153F5.1) 
WRITE( 6310) MsNoR 
Mi=M+1 
Ni=N+1 
MM1=M—1 
NM1=N-1 
VZB=1.0 
READ(5.15) DsRC 
FORMAT (2F10.0) 
READ(5s16) ITMsISMsIVM 
FORMAT( 3110) 
READ(5317) OMEV,sOQMESsO0MEW 
FORMAT(S3SFi0.1) 
READ(5,117) FRI1 
FORMAT(F10.3) 
ESS=1.E-4 
PSS=1.E-3 
HI=(1.4+1./R)/7( 20M) 
HJ=(1.#+R)/S(4.*N) 
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HI2=HI*HI 
HJ2=HJ*HJ 


HO=-2 .*(1./HI24+1./HI2Z) 


DO S I=1,Mi 


XCT)=-01.241./R2) 44.40 1-1) *HI 


RA1(1)=1.4+X(1I)/RC 
RA2 (I )=RA1 (1) *RAIC 1) 
CONTINUE 

DO 450 I=1,.Mi 


READ(5s460) (SC IsJ) sJ=19N ) 


DO 451 I1=1.M1 


READ(53s460) (WCIsJ) sJ=19N ) 


DO 452 I=1>M1 


READ(5s460) (VOR( 15S) sJ=15N1) 


DO 888 I=1.Ml1 


READ(52460) (UCIsJ) »J=1.N1) 


DO 889 I1=1,.,Ml1 


READ(5s460) (VC IsJ) sJ=19Ni1) 


FORMAT (5E16.7) 
BOUNDARY CONDITUONS 
DG 110 I=1.Mi 
S(I.1)=0. 

SC IsN1)=0. 
VCIst )=0¢6 

VC IoN1)=0. 
UCIsN1)=0. 
wCIsN1)=0. 

DO 111 J=i,N1 
SOils J3=06. 
UGiss}=0. 
U(M1sJ)=C. 
V(M1sJ)=0. 
VOLO I VEC < 
S(M1.J)=0. 
WCM1,J)=0. 
WOERSOUDSEC. 
VGR(1-eJ)=0. 
Ga 

CONTINUE 

DO 12 I=2>5M 
DO 12 J=2N 


UC Is J) =0.5*(SCIsJt1)-SC155-1)) 7 CHI *RAI CI) ) 
VOIsJ)H065*(S(T-1 J) -SCIt1s J) IS CHI*RAI CI) ) 
DW(IsJ)=W( I oJ) *(WC Is Jt1)-W0 1s 5-1) 7 (RC*XRAZC(I) *HJ) 


DO 13 I=25M 


UC1Ts1)=0-5¥*(4.*S(152)-S(153))/(HI*RAI(I)) 
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VORTICITY AT BOUNDARIESP 
DO 20 J=1.Ni 
VOR (1 oJ) =—-26*S (295) ¥VZB/(RA2(1) *¥HI12)4(1.-VZ8) ¥VOR(12J) 
20 VOR(M15IS)=-2. *S(Ms J) *XVZB/(RA2Z(M1) ¥H12) 4(01.2-VZB) ¥VOR(M15J) 
DO 21 F=1.M1 
21 VORCIsN1)=-2.*S( IN) *¥VZB/(RAZ( I) *HJ2)4+(1.-VZ8) * 
iVOR(IsN1) . 
SOLVE VORTICITY FUNCTION 
Iv=1 7 
160 CONTINUE 
DO 22 I1=2.5M 
DO 22 J=2.N 
22 S1(IsJ)=VOR(IsJ) 
DO 25 J=25N 
DG 23 I1=2s5M 
L=I-1 
A1 (L) =€1./HI 2-1 -5/( RCXHIXRAI (I) )F0.5*UCToJ)/SHI) 
A2(L)=HO 
AZ(L)=(1-/HI 241. 5/(RCXHI*¥RAI1(1))-0.5*UCI5J)/HI) 
A4(L)=DW( I oJ) —-(1./HI2-0 oS *¥VC Is J )/SHI)D XVOR( I, 541) 
* —-(1 SHI2460 -S*V(C Is JS) /HI) *VOR(CI,J5-1) 
23. CONTINUE 
A4 (1) =A4(01)-A1(01)*VOR(1sJ) 
A4 (MM1)=A4 (MM1 )-A3(MM1 )*VGR(M1 oJ) 
CALL GAUSS(A1.A2sA33A45MM1 5X) 
DO 24 I1=2.5M 
VX(I5J)=xX(iI-1) 
24 VOR(1I.J)=S1I(C 16S) tOMEV*X(VX(I 5S) —-SIC IJ) ) 
25 CONTINUE 
DV=O4 
FEV=0. 
DO 30 I=2M 
DO 30 J=2:N 
VN=ABS(VOR(I6J)) 
VR=ABS(VOR(I5J)-S1(IsJ)) 
DV=AMAX1(D0V.2VR) 
FE V=AMAX1 (FEV>sVN) 
30 CONTINUE 
DV=DV/FEV 
IF (DVeLE-ESS-e-OReIVeGT.IVM) GO TG 31 
760 Iv=IVvtli 
GO TO 160 
31 CONTINUE 
SOLVE STREAM FUNCTION 
DO 310 I=25M 
DO 310 J=2sN 
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300 


32 


34 
35 


40 


41 


42 
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S2(1,J)=S(1 oJ) 

ISs=1 

CONTINUE 

DG 32 I=2>5M 

DO 32 J=2sN 

S1(IsJ)=S( 15d) 
DQ 35 J=25N 

DO 33 I=25M 

AVP=-RA2(I)¥*VOR(I4J) 
L=I-1 a: 

A1 (LJ) =(1e/ HI 240 65/ (RCXHI*RAI(1))) 
A2(L)=HO 

AZ (L)=(1./HI 2-0. 5/ (CHI*RC¥RAI(I))) 

AG (L)=AVP—(1 -/HI2Z) (SC IsJt1)4+S(15,5-1)) 
CONTINUE 
CALL GAUSS(A1+sA2sA3 5A49MM1.X) 

DO 34 I=25M 

VXCIsJS)=XCI-1) 
S(IsJ)=S1C1sJ)4+OCMES*(VX(IsJ)-S1(1I9J)) 
CONTINUE 

HS=oe 
FES=0. 

DO 39 I=2.5M 

DO 39 J=25N 

SN=ABS(S(1I2J)) 

SR=ABS(S(CI9J)-S1I(IsJ)) 
DS=AMAX1(DSSR) 

FES=AMAX1(FES>3SN) 

CONTINUE 

DS=DS/FES 

IF (DS.LE.ESS.~OR-IS-eGTLISM) GO TO 40 
IS=ISti 
GO TO 300 

CONTINUE 

WRITE(65,350) IS,DS 

UP=-DATING COEFF 

DOSS MALS 2 9M 

DO 41 J=25N 

UC Ts J) =0.5*(S(1 5Jt1)-SC Is J-1)) /7(HJ*RAI (I) ) 
VOITeJ)=0.5*(SC I-19 J) -SCT#t1sJ)) /CHI*RAICI) ) 
DO 42 I1=25M 
UCTs1)=0.5*(4.*S(1s2)-SC193))/(HI*RAICI)) 
DSs=0). 

DO 711 I=1.M1 

OG EF EY HISRSN1 

SRR=ABS(S( 155) -S2C1595)) 
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5S 


5000 
54 
55 


59 


64 


350 


DSS=AMAX1(DSS+sSRR) 
CONTINUE 
SD=DSS/FES 
SOLVE AXIAL VELOCITY 
DOW71L2AT=25M 
DON? IT 20g=2.N 
S2(IsJ)=Ww(I oJ) 
Iw=1 
CONTINUE 
DO 51 1=2.5M 
DO 51 J=25N 
S1(TsJ)=W(IsJ) 
DG S2 I=25M 
WCI15,1)=(4.*wW(1,2)-w(1,3))73. 
DO 55 J=2,N 
DO 53 I=2.M 
L=I-1 
Al (LJ) =€1 «-/HI 240 .5%*UC I eo JJSHI-O65/(RCXHI*¥RAI(I)) ) 
A2 (L)=(HO-1 -O*U(IsJ)/0RCKRAI (1) )-1.0/(RC*ROXRA2Z(I))) 
AZ (LJ=(12/HI2-0.5*UC Is JJSHI4F0 -5/ (HI *RC¥RAICI)) ) 
C=D/RAI(I) 
A4 (LL) =C—(1.6/HI 2-0 SVC Ie J) SHI) KWOT S Jt1) -C 1. SHI 2+ 


10.5*V( IJ) /HI)*WC15J-1) 


CONTINUE 

CALL GAUSS(A1 sA25A3 sA4 9 MM1 2 X) 
DG S4 I=25M 

VX(IsJ)=XCI-1) 
FORMAT(11E10.3) 
W(IleJ)=S1(0 155) tOMEWX(VX(I 65) -Si(I2J)) 
CONTINUE 

Pw=0. 

FEW=0. 

DO 59 I=25M 

DG 59 J=2sN 

WN=ABS(WC(12J)) 
WR=ABS(W(12sJ)-S1(IsJ)) 
PW=AMAX1(WRPW) 
FEW=AMAX1(FEWsWN) 

CONTINUE 

PwW=PW/FEW 

IF (PW.LE.ESS.~OGReIWeGT.ISM) GG TO 64 
Iw=Iwtl 

GO TO 750 

CONTINUE 

WRITE(65350) IW+PW 

FORMAT (/s153E20.5) 
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Pww=0. 
DO 714 I=26M 
DO 714 J=2s5N 
WD=ABS(W(1IsJ)-S2(1sJ)) 
PWW=AMAX1 (WD »Pww) 
CONTINUE 
WO=PWW/FEW 
WRITE(63;223) SD.wWD 
FORMAT(*SD =",615.5s/e"WD =',E615.5) 
WRITE(6s222) K 
FORMAT(*OQUTER ITERATION TIMES =',415) 
FORMAT (15915XsE15259//915915XsE15.5) 
DO 65 1=25M 
DO 65 J=2.N 
DW(IsJJ=WC 1 oJ) *CWC I sJt1)-WC 1s J-1))/ CRCXRA2(T) XH) 
IF (WD.LEePSS.AND.SD.LE-PSS) GO TO 200 
IF(K.GEsITM) GO TO 200 
K=K+1 
GO TO 660 
CONTINUE 
CONTINUE 
CALCULATION OF DEAN NUMBER 
WRITE (6310C5) 
FORMAT(® STREAM FUNCTION®) 
DO 1006 I=1.M1 
WRITE (65601) (SC IoJ) »J=15N) 
DO 3001 I=1.M1 
WRITE(6-601) (VOR( I» J) eJ=15N1) 
CONTINUE 
WRITE(6 61) 
FORMATC®AXIAL VELCCITY®) 
DO 600 I=1.M1 
WRITE(62601) (WC IsJ)sJ=15N) 
WRITE(691001) 
FORMAT (C®V COMPGNENT ® ) 
OG 1177 T1l=1M1 
WRITE(6.6C1) (VC IsJ) »J=1 oN) 
WRITE( 651002) 
FORMAT(*®*U CGMPGNENT ® ) 
OGIt73 s1=1 71 
WRITE(6s601) (UCI sJ) sJ=12N) 
CALL SIMPS(WoMesNsHI sHJs SUMM) 
AREA= (1.4+R)*(1.4+1.6/R)/8 
SUMM=SUMM/ AREA 
WRITE(62501) SUMM 
FORMAT(*REYNOLD NUMBER =" 510X5,E15.5) 
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SUPP=SUMM/SQRTC(RC) 

WRITE(631000) SUPP 

FORMAT(//+s"*DEAN NUMBER =", F20.5) 

DO 130 I=1.5Mi1 

FIC I) =-0.5%*(-4.*W( I aN) t+WOITSNM1))/7 HJ *RAI(I) 
CALL SIMP(CFJsMsHI.SUM) 

SUM1=SUM 

DO 132 J=1,N1 

FIC JI )=-06.5*(0-4e*¥ WC Ms J) EWC MM1 55) ) SHI *(120¢+0.2.5/RC) 
CALL SIMP(CFJsNoHI SUM) 

SUM2=SUM 

DO 135 J=1,N1 

FJ (5) =-0.5%(-4.*W( 205) tW0 Sed) ) SHI *¥(12-9657RC) 
CALL SIMP(FJsNoHIs,SUM) 

SUM3=SUM 

PRE=0.5*(1.+R)4+0.5*(1.4+1./R) 
FRE=(SUM14+SUM24SUM3 )/PRE 

FRE1=2.*FRE/SUMM 

FRAI=FRE1/FR1 

FRE2=-0.5*D/SUMM 

FRA2=FRE2/FRI1 

WRITE(6s666) FRE1,FRE2, FRA1»s FRA2 

FORMAT(CS// so *FREL ="6F15.595Xe*FRE2 ="'sF1S.59// 
*K*®FRAL =*,F152.595Xs"FRA2 =",F15.5) 

FORMAT(S »+7E18.5) 

FORMAT(*® *,13310E12.3) 

DOl4ZO151=9 SMi 

WRITE(72410) (SC IsJ)sJ=19N ) 

pO F402 TI si TMs 

WRITE(72410) (WC ITesJ)sJ=15N ) 

DO 403 I=1.Mi 

WRITE(75410) (VOR( Is J) »J=15N1) 

DG 190 I=1.M1 

WRITE(72281) (UC IsJ)sJ=15N1 ) 

DO 191 I=1.M1 

WRITE(72281) (VC IsJd)sJ=19N1) 

FORMAT (5E16.7) 

FORMAT(5£E16.7) 

CALL CCC2CITIsJSsMAXTsXTsTsHIsHJ) 
WRITE(62531) MAXTsIIsJJ 

CALL CCC1CXTsTsHIsHJ) 

FORMAT (*0%s20Xe*MAXIMUM OF TEMPERATURE ® 910K s£14.7 
195Xs*IN I-D*s2XeI5e5Xe"IN J-D's2X5I5) 

STOP 
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SUBROUTINE FOR THE SELECTING GF THE MAXIMUM VALUE 
1 IN THE NUMERICAL CALCULATION DOMAIN 
SUBROUTINE CCC2¢(11sJJs MAXTSsDEsSTsHIsHy) 
DIMENSION ST(41521),DE(10) 

COMMON /C3/ MN 

REAL MAXTS»sMAXST 

FES=0. 

DO 1 I=2.5M 

DO T19I3=23'N 

SC=ABS(ST(1I4J)) 

FES=AMAX1(FES3SC) 

DO 2 I=2.5M 

DO 2 J=25N 

SC=ABS(ST(1I5J)) 

IF (SC.EQ.FES) II=I 

IF(SC.CEQ.FES) JJ=J 

CONTINUE 

MAXTS=FES_ 

MAXST=FES 

EH=0. 

NJ=1 

IM=2 

IF (MAXST.LT~20.00001 .OR.MAXST.GT.2100) GO TG 50 
IF (MAXST.GE.0200001 «-ANDeMAXST2LT20.0001) DB=0.00001 
IF (MAXST.GE.0.20001 «-ANDeMAXST.LT202001) DB=0.0001 
IF (MAXST.2GE.20.001 eANDSMAXST.LT.0201) DB=9.001 
IF (MAXST.*GE.0.01 eANDSMAXST.LT~20.01) DB=0.01 

IF (MAXST.GE.061 eANDeMAXST.LT~21.0) DB=0.1 

IF (MAXST.GEe1 -0 eANDSMAXST.~LTL10.) DB=1.0 

IF (MAXST. GE.10 «eAND.MAXST.LE~2100.) OB=10. 
DA=DB 

DO 4 1=2s9 

DM=AMIN1I(MAXST DA) 

IF (DMZEQeDA) DD=DM 

DA=DBx I 

CONTINUE 

EH=OD+EH 

IF(NJ.2EQ.2) GG TO 8 

MAXST=MAXST-DD 

NJ=NJ4+1 

GO TO 7 

DC=EH/10. 

DONS) J=1.10 

DE (J) =DC*J 

CONTINUE 
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IM=1 
RETURN 
END 


SUBROUTINE FOR THE CONSTANT STREAM LINE 
SUBROUTINE CCC1 (CFeFsHIsHJ) 
DIMENSION CF(10).F(41.21)5D0F(41)+.CI(41)5PF(41 ) 
2R1I1(41).PI(41) > Q1I(41) 
COMMON /C3/ MoN 

DH T62K=1.10 

WRITE (69120) KsCFC(K) 
WRITE(65105) 

DO 81 I=25M 

JK=1 

feEhSGE S2)) “GO *TOe>Ss 

DF (1) =CFCK)-F(CI51) 

DG 83 J=2.,N 

DF (J) =CFCK)-FCI.J) 
EE=D0F (J) *DF(J-1) 

IF(EE.GT.O.) GO TO 83 

EF=DF (J-1)/(F (I1sJ)-FCIsJ-1)) 
CJC JK )=J—-1+4EF 

JK=JK+1 

CONTINUE 

IFC JK.EQ.1) GO TO 81 

JKD=JK-1 

QI (1) =HI*CI-1) 

WRITE(62107) IsQI(1I)s(CJC JK) sJK=15d5KD ) 
CONTINUE 

WRITE(62108) 

DO 71 J=2sN 

i= 1 

IF(J.GEe-2) GO TO 74 

PF (1)=CFCK)—-FCisJ) 

DG 72 I=25M 

PF (I )=CF(K)—-FCIsJ) 
QF=PF(1)*PFC(I-1) 

IF(QGFeGT.O.) GO TO 72 
RF=PF(I-1)/7°0FCIeJS)-FCI-15J)) 
RIC IK )=1-1+RF 

IK=IK+1 

CONTINUE 

IF CIK.EG.1) GO TO 71 

IKD=IK-1 
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PJ(J)=HU*(J-1) 
WRITE(691607) JsPJICJ)oCRICIK) » IK=1.I1KD) 
71 CONTINUE 
6 CONTINUE 
105 FORMAT( 80%, 30Xe*COGRDINATE IN X—-DIRECTIGON' »20Xo°!I 
IN Y-DIRECTION®) 
107 FORMATC® ® ,20X%s1595X9F1004920X%54F10.4) 
108 FORMAT(#0*,30Xs*CGORDINATE IN Y-DIRECTIGON's20X>o 
1 *X-DIRCETIGN® ) 
120 FORMATC #0" »s20Xe"*NO GF CONSTANT LINE" s10X515010X» 
-1 "CONSTANT VALUE® 510XsF20.5) 
RETURN 
END 


KK ok ak SK ok ok ok sk ok ok ek ok ok SKK ok 3K Be ok SK Sk Ke OK OK 3K 3K KK RK Ok OK KE OK ak OK ie OK OK KK KK KKK 
HEAT TRANSFER OF UNIFORM WALL TEMPERATURE 
3K XK ok SK ok ok oe ok ok 3 ok ok He Ke ek KK KKK a a SK He ae eK RK KK KK Ke eK IK KK KK KK OK A OK KK 
R ASPECT RATIG 
REY REYNOLD NUMBER 
DEAN DEAN NUMBER 
PR PRANDTL NUMBER 
DIMENSION 17(41421) 571041521) 5TK(41521) »7T2(041521) 
DIMENSION A1(4i1 521) .A2(41.21 )+A3(41521) 5B1(41921)s5 
*B2(41521)+.63(415,21) 
DIMENSION $(41521)sVGR(41.221) »sFJ(41) 
DIMENSION X(41) »sRA1(41),SWT(2090) 
DIMENSION U041521)5V(41521)5W(41521) 
DIMENSION GZW(41 521)eWT(41521) »S0Z(200) 
DIMENSION XT(41 21) 
COMMON /C1/ ZZ(1500)sDXsKZ » BI 
COMMON /C2/ Al4121)sBsC( 41521) 20(41521)3EsF( 41.21) 
1ww(41.21): 
COMMON /C4/ M1;5N1 
COMMON /C3/ MsN 
COMMON /C6/ AREA 
REAL MAXS »MAXT 
READ(Ss10) MsNoR 
WRITE(6210) MsNoR 
10 FORMAT(2155F5.1) 
READ( 5s 314) PRsRC»DEANSREY 
314 FORMAT(4F20.5) 
WRITE(6014) PRsRCsDEANSREY 
14 FORMAT(//s*PRANDTL NUMBER ='9F15.59/s5 
1*CURVATURE ="sF15.59/ +s 
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**DEAN NUMBER =* 6F15.59/s"REYNOLD NUMBER =*5F15.5) 
Mi=M+1 
N1i=N+1 
MM1=M—1 
NM1=N-1 
HI=(1.t1./R)/(2.%*M) 
HJ=(124+R)/S(4.N) 
DX=HI/2. 
DX 2=HI*HI 
DO 450 I=1.M1 
450 READ(5s460) (S(IsJ) sJ=1oN ) 
DO 451 I=1.M1 
451 READ(55460) (WC IsJ) sJ=15N ) 
DO 452 I=1>5Mi1 
452 READ(5:460) (VOR(1I5J).s J=1 »N1) 
DO 888 I=1.M1 
888 READ(5s460) (CUCT.J) sJ=15N1) 
DO 889 I=1.5Mi 
889 READ(5+460) (VC IsJ) sJ=15N1) 
460 FORMAT(5E16.7) 
ENTRANCE CONDITIONS 
00 £5 1=1.1 
w(IseNi)=0. 
DO 15 J=1.Ni 
15 TCLSY}=150 
DO 16 I=1,M1 
16 TCI+sN1)=0. 
DeS26390=15N1 
TCLs SI=Gs 
26 ©TCM1 PT I)H08 
DO 12 I=1.Mi1 
H0012 359 +15N1 
12 WCIsJ)=WCIsJ)/REY 
DO 100 I=1,.Mi 
XCT)=—- C1 .t16/R)/4e4( 1-1) XHI 
100 RA1(1)=1.t4+XC1I)/RC 
DO 18 I=15M1 
DO 18 J=1,N1 
18 WWC(Is/J)=2.*W(IsJ)*DX2/RA1(1I) /PR 
PR1=1./PR 
DO1t~ 1T=15M1 
DG 11 J=1,N1 
AC(IoJ)=PR1—-DX* (UCT sJ)—-1./(RC*XRAICI)) ) 
B=-2.*PR1 
C(le J )=PR1I4tDX*(UCTsJS)—-1./(RC#RAICI)) ) 
D(C IsJ)=—-PR14+DX¥*¥V(I5J) 
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E=2.*PR1 
FCI eJ)=—-PR1I-DX*V(I5J) 
CONT INUE 
AREA=(1.4R)*(1.4+126/R)/8. 
KMAX=1000 
Z=0. 
KZ=0 
KK=1 
BOUTS =T=PIM1 
DO 19 J=1,.N1 
TK(CIsJ)J=TC1IsJ) 
DX=HI 
KZ=KZ+1 
AXIAL STEP 
DZ=5.E-5 
IF (KZ.GT.2100) OZ=1.E-4 
TF (KZ.GT.450) DZ=2.E-4 
IF(KZ.GT.650) DZ=5.E-5 
IF (KZ.GT.850) DOZ=1.E-3 
Z=Z+DZ 
ZZ(KZ)=DZ 
DO 22 I=1.M1 
DO 22 J=1,5N1 
TKCIsJ)=T(I5J) 
CALL TEMP(T »TK) 
WRITE(62201) KZsZ 
DO 23 I=1.M1 
DO 23 J=1;N1 
TICIsJ)=H=TC I oS) *¥WOIT Ss J) 
CALL SIMPS(T1s+MsNsHIsHJ,SUMM) 
TBK=SUMM 
DO 136 I=15Mi 
FUCT)=-0.5%*(-4.*T( IsN)+TCI2NM1))7 HJ 
CALL SIMP(FJsMsHI»SUM) 
SUM1=SUM 
DOT 1S 21 9=1 sN1 
FIC J) =-0 6 5*(-4 0*T (Ms 5S) 4+T(MM12J5)) SHI 
CALL SIMP (FJUsNeHJ»SUM) 
SUM2=SUM 
OO £1335 I=15N1 
FJ (J) =-0 65 *(-40.*T( 295) 470395) ) SHI 
CALL SIMP (FJIsNoHJsSUM) 
SUM3=SUM 
PRE=0 .5*(1.64+R)40.5%(1.4+1./R) 
SUMM= (SUM1+4+SUM24+SUM3 )/PRE 
ANU1=—-SUMM/T BK 
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CALCULATIONS GF TEMP. GRADIEND IN FLOW DIRECTION 


DO 5O I=1.M1 
DO SO J=1.Ni1 
T1CIsJJ=H( TCI 9S) H-TKCIs J) )/SZZ(KZ) 
DG 51 J=1,N1 
DO 51 I=1.Mi 
TICIsJ)=H=TICIs J) *WCI oJ) 
CALL SIMPS (T1 o«MsNeoHIsHJs SUMM) 
ANU2=-0.25 * SUMM/TBK 
AVNU=0.5%*( ANUI1+ANU2 ) 
DNU=100*ABS (ANU1—ANU2) /AVNU 
IF(KZ.LE.1) GO TS 66 
BNU2=0.5*( ANU2+BNU2 ) 
BNU =0.5*(BNU14B8NU2) 
BP=100. *ABS (BNU2-BNU1 )/BNU 
WRITE(6202) BNU1+s8NU2 «BNU, BP 
CONTINUE 
BNU1=ANU1, 
BNU2Z=ANU2 
WRITE( 62501) ANU1sANU2sAVNUSDNU 
WRITE(6,502) TBK 
WRITE(69503) 
FORMAT(/s*NG. OF SECTIGN's15310Xs "DISTANCE FROM 
1ENTRANCE® .F10.6) 
FORMAT (30X 5 **** ",3F15.63£20.5) 
FORMAT ( *"NU1 * 6F10.6310X9 "NUZ2%sFi0.55910X%s "AVERAGE 
1NU "5F10.6510XsE2004) 
KKK=KZ/30 
IF (KKK.NE.KK) GO TO S00 
KK=KK+#1 
WRITE(63503) 
DO 61 I=1.M1 
WRITE(65504) I 
WRITE(65505) (TC IsJ)sJ=12N1) 
FORMAT(*®* BULK TEMP.'»s £15.25 ) 
FORMAT(/S/ *TEMP. DISTRIBUTIGNS!/) 
FORMAT(*I*,15) 
FORMAT( * *,10612.4) 
CONTINUE 
IF (KZeGEeKMAX) GO TO 30 
GO TO 20 
CONTINUE 
STOP 
END 
SUBROUTINE TEMP(TsTK) 
DIMENSION T(41 521) sTK( 41521 )sTKH(41521) 5X(41) 
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DIMENSION A1(41) .A2(41 )3A3(041)5A4(41) 5A5( 41) 
COMMON /C1/ 2Z(1500) sDXsKZ » BI 
COMMON /C2/ A(41921)sBsC(41521) 5D0(41521)5EsF(41521) 
1.ww(41.21) 
COMMON /C3/ MsN 
COMMON /C4/ Mi1,5N1 
M1i=M+1 
MM1=M-1 
DO 20 J=1.N 
DO 10 I1=2.5M 
L=I-1 
A1(LJ=CCI sy) 
A2(L)=B-wWW(IsJ)/ZZ(KZ) 
A3(L)=A(1I.J) 
fe CIVEQS 1) “GOTO TS 
A4(LJ=D(C Io JS) *XTK(C Io S41) HC E-WWO ITs J) SZZ(KZ)) *¥TK(C Iss) + 
LFCIsJ)*TK(CIsJ5-1) 
GO TO 10 
A4(L)J=( OCI oJ) FFCI 5S) )*X¥TKC Is S41) 40 E-WWC led) ZZ(KZ)) 
1*TK(I».J) 
CONTINUE 
CALL GAUSS(A1e¢A2sA3sA49MM1 2X) 
DO 20 I=25M 
TKHC Is J)=X(I-1) 
DG 21 J=1.N 
TKH(1.2J)=0. 
TKH(M15J)=0. 
DG 40 I1=2.5M 
DO 30 J=1+N 
Al CSV=F(1-.J) 
A2(CI)=EtwwlsJ)/7ZZ(KZ) 
A3(J)=D( 10d) 
AG (J) =ACT os J) XTKHC I 415/39) 4° BtWww( lI sJ)/ZZ0KZ) ) *T 
1KH( 16d) 4CCI sd) XTKHC I-15 J) 
CONTINUE 
ASCL)=FCIs1V+D (Ist) 
CALL GAUSS(A1 sA2sA3 2A4 No X) 
DG 40 J=1.5N 
TCIsJ)=XCJ) 
RETURN 
END 
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Sak ook aK ak ok ok aK ok kok 2 ook fei ok ake ok ak ake ok Sa akc oi ako kk ak ok ak ak ae a fe 2 ofc ac ok fe ak ok akc ak ok 3k ak akc aie ak ok 
R ASPECT RATIO 
REY REYNOLD NUMBER 
DEAN DEAN NUMBER 
PR PRANDTL NUMBER 
DIMENSION 1T(41 421) 971(041521)5TK(41.21) »T2(41521) 
DIMENSION S$(41421)sVOR(41521) +FJ(41) 
DIMENSIGN X(41) 5RA1(041).SWT( 200) 
DIMENSION U(41521)5V(41521 )sW(41,21) 
DIMENSION GZW(41 621) »WT(41521) »sSO0Z(200) 
DIMENSION XT(41 521) 
CGMMON /C1/ ZZ(1500).D0X+KZ 
COMMON /C2/ Al 41521) 1€041,21):0(41,21) 3F(41,21) 
oWW(41 4.21) .8B(41521) 5E( 41,21) 
COMMON /C4/ MiNi 
COMMON /C3/ MeN 
COMMGN /C6/ AREA 
REAL MAXS»sMAXT 
READ(5310) MsNoR 
WRITE(6910) MsNoR 
FORMAT (2156F5e1) 
READ(5s314) PR»RCsDEANsREY 
FORMAT(4F20.5) 
WRITE(6514) PR»eRC»sDEANs REY 
FORMAT(//>"*PRANOTL NUMBER =",5F15.50/5 


1*CURVATURE ="sF15.59/s 
*"*DEAN NUMBER =" 6F15.59/7s"REYNOLD NUMBER =! 5F15.5) 


Mi=M+1 

Ni=N+1 

MM1=M—-1 

NM 1=N-1 

HI=(1.+1./R)/(2-*M) 
HJ=(1.24+R)/(4.N) 

DX=HI/2. 

DX 2=HI*HI 

ENTRANCE CONDITIONS 

DO 15 I=1>5Mi 

WC IsN1)=0. 

DE4tS J=LI9NA 

TCI sJ)=0.0 

DX=HI 

DO 450 I=1,5Mi 

READ(5,460) (S(IsJ) sJ=15N ) 
DO 451 I=1.M1 
READ(S5s460) (WC IeJ) sJ=15N ) 
DO 452 I=1.M1 
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READ(52460) (VOR(CIsJ)s J=1 5N1) 
DG 888 I=1.M1 
READ(5s460) (U(IsJ) sJ=15N1) 
DO 889 I=1,.M1 
READ(53460) (VC IsJ) »J=12N1) 
FORMAT(5E16.7) 
DO 100 I=1,M1 
X(T )=-€1 2.41 ./R)/4.4+0 1-1) *HI 
RAI (1 )=1.4+X(1)/7RC 
DO 12 I=1.Mi 
DP 12) JHE a NI 
WCITsJ)=WCIsJ)/REY 
DO 18 I=1.M1 
DO 18 J=1,.N1. 
WWCIsJ)=2.0*W( I 9d5)*DX2/RAI (I) PR 
PR1=1./PR 
DO 11 I=1.5M1 
DGO411 JJ=1 sNI 
AC IsJ)=PR1-DX*(UC 1, I)-16/(RC#RAI(I1)) 
BCIsJ)=-2.*PR1 
CCI sJ)=PR140X*X(UCTsJ)—-1./7(RC#¥RAI(I)) 
D(C IsJ)=—-PR1+DX*V( IJ) 
E(I.,J)=2.*PR1 
F(IsJ)=-PR1-DX*V( IJ) 
CONTINUE 
AREA=(1.4R)*( 1.241 .6/RI/86 
KMAX=1000 
Z=0 6 
KZ=0 
KK=1 
KZ=KZ+1 
AXIAL STEP 
DO 19 I=1>5M1 
DO 19 J=1.N1 
TK(Is J) =TCI6J) 
AXIAL STEP 
DZ=5.E-5 
IF (KZ~.GT.100) DZ=1.E-4 
IF CKZ3GT3450) (DZ=2E=4 
IF (KZ.GT.2650) DZ=5.E-4 
IF (KZ.GT.850) DZ=1.E-2 
Z=Z+DZ 
ZZ (KZ) =0Z 
DO 22 I=1.M1 
DO P22 .J=15N1 
TK(C Te J) =HTCI SJ) 
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CALL TEMP (T.TK) 
WRITE(65,201) KZsZ 

DO 23 I=1.M1 

DO 23 J=1.N1 
T1ICIsJ)=TCI sJ) *XWCI oJ) 

CALL SIMPS(T1 »sMsNsHIsHJsSUMM) 
TBK=SUMM 

SUT=0. 

DO 42 I=15Mi 

SUT=SUT#TCIN1) 

DO 44 J=1eN 
SUT=SUTtT(1sJ5)4+T(M1 oJ) 
SUT=2.* SUT-T(151)-T(M1 o1) 
AVET =SUT/C2*(M424N) ) 
TAW=AVET 

DO 651 I=1.Ml1 

DO 651 J=1;5N1 

TICIT se J) =(AVET—-TCIsJ))*¥WCI J) 
CALL SIMPS(T1 sMeNvsHIsHJ»SUMM) 
ANU1=1./SUMM 

TMW=SUMM 


CALCULATICNS OF TEMP. GRADIEND IN FLOW DIRECTIGN 


DO 25 J=i>sN1 

DO 25 I=1,.M1 

GZW( Ile J) =( TCI 5s SJ) -TKCI55)) *WOI5 5) /7ZZ(KZ) 

CONTINUE 

CALL SIMPS(GZW »sMsNeHIsHJs SGWT) 

ANU2 =0.25 * SGWT/TMW 

AVNU=0.5%* (ANU1+ANU2) 

DNU=100.* ABS(ANU1—-ANU2)/AVNU *0.25 

WRITE( 65501) ANU1 » ANU2 » AVNUsDNU 

WRITE(6.502) TBK sTAW 5TMW 

KKK=KZ/30 

IF (KKKeNESKK) GO TO 500 

KK=KK+1 

WRITE(62503) 

DO 61 I=1.M1 

WRITE(65504) I 

WRITE (62505) (TCIsJ)sJ=15N1 ) 

FORMAT(/s"*NGe OF SECTION! »15910Xs*DISTANCE FROM 
1 ENTRANCE® »F1i0.6) 

FORMAT (30X%»9 *****,3F15.63£20.5) 

FORMAT( ®NU1 *§ 5F10.6910Xs *NU2Z "9 F10.6910Xs "AVERAGE 
1NU "sF10.66310X3E20.4) 

FORMAT(*BULK TEMP .%s £15.65 s*WALL TEMP.'sE15.5s 
1"MEAN WALL DIFFERENCE*, £15.25) 
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FORMAT(// "TEMP. DISTRIBUTIGNS'Z) 
FORMATC*I*®.1I5) 

FORMAT( * *,9£12.4) 

CONTINUE 

IF (KZ~GEeKMAX) GO TO 930 

GO TO 20 

CONTINUE 

SToP 

END 

SUBROUTINE TEMP(T.TK) 

DIMENSION 1T(41 921) sTK( 41921 )sTKH(41521) »X(41) 
DIMENSION A1(41)sA2(41)sA3(41) 5 A4(41) 5A5(41) 
COMMON /C1/ ZZ(1500)s»DXsKZ 

COMMON /C2/ A(41s21) 31€(41421).D(41,21) »F(41,21) 


16WW(41521) .8(041221).E( 41.21) 


COMMON /C3/ MsN 
COMMON /C4/ Mi1sN1 
DO 20 J=1;3N 
DD 10 I=2—M 
Al (1)=CC1sJ) 
A2(I)=BCI1sJ)-WWO Ts J) /ZZ(KZ) 
A3(I)=ACI +4) 
IF (JeEQe1) GO TO 15 
AG(I)=D( Is J) *TKC Is Jt1) +0 EC I, 5) -Ww(leJ)/ZZ(KZ) )*TKCI 9d) 


*4F (I1.J) *TK(IsJ-1) 


— 


GG TO 10 
AG (I)=(OCI sJPtFFC I oJ) P¥TKC Io JFt1L FC ECIs JS) -WW(I 5s) 
JZZ(KZ))V*TKCIsd) 
CONTINUE 
A3(1)=2. 
Al Gat }=2. 
A2(M1)=-2. 
TE CESEQSE)LMGD TO 16 
A4(1)=—TKC 15541 )426.*TK( 19S) —-TK( 155-1) -—2.*DX 
AG OM1)=—-TK(M15J41)4+2.*TK( M195) -TK(M1sJ-1) -2.*DX 
A401) =-2.*TK(1 52) 42.*TK(151) -2.*DX 
A4 (M1) =—-26*TKCM1 92) 426 *TK (M191 )—-2- *DX 
CALL GAUSS(A1lsA2sA3sA45M1 9X) 
OGEY7en—15M1 
TKH(C Is J)=X(I) 
CONTINUE 
DO 21 I=1,.Mi 
A1(I)=1. 
A3(1I)=1. 
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A4(1)=-2.*TKH(IsN)—-2.*DX 
A3(1)=2. . 

A1 (M1 )=2. 

A4(M1)=A4(M1)-2.*DX 
A4(1)=A4(1)-2.*DX 
CALL GAUSS (A1sA2sA3sA49MizsX) 
DO 22 I=1.M1 
TKHCIsN1)=X(T) 
DO 40 I=2.5M 

DO) 30 -J=1,5N 

A1(J)=FCIsJ) 

A2C J) =ECIs J) tWWOIT SI) /ZZ(KZ) 
A3Z(J)=DC1I5I5) 

A4(JSJ=ACT so JD ¥TKH(C 14154) 40 BC Ilo J) tWWCIs JIS ZZ(KZ) XT 


*KH(CIsJ)4+CC 1,5) *TKHCI-15J) 


CONTINUE 

AITON1)=2. 
A3(1)=F(1I151)4+D(1I51) 
A4(N1)=—-TKH(CI+19N1 )4+2.*TKHCIsN1 )-TKH(CI-15Ni)—-2.*DX 
CALL GAUSS(A1sA2sAZ3 sA49N1 9X) 
DO 40 J=1.N1 

TC Is J) =XCJ) 

DG SO J=1,Ni 

Al(J)=1. 

A2(J)=-4. 

ASCS)P=1 6 


(AG (I) =—-26*T(295)—-26 *DX 


AS (J) =-2.*T(MsJ)—2.%*DX 
A3(1)=2. 
A1(N1)=2. 
A 4(N1)=A4(N1 )-2.*DX 
CALL GAUSS (A1sA2sA3sA45N1 9X) 
ba 'S2+3=1.N1 
TCLs J) =X( J) 
A5 (N1 )=AS5CN1 )-2.*DX 
CALL GAUSS(A1sA2sA3sA59N1 5 X) 
DO (54 \J=1.N1 
T(M1sJ)=X(J) 
RETURN 
END 
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SKK He Sk aie SKK OK Oe ok SK Sik aie OK OK ok a ok sk Se ak 3c ok AK KK OK aK ok Sk aK OK a 3 Kk i kK ok ok OK KKK KK KKK KK KK 
R ASPECT RATIO 

REY REYNOLD NUMBER 

DEAN DEAN NUMBER 

PR PRANDTL NUMBER 

BI BIOT NUMBER 

DIMENSION 1(414521)2T1(041.21),TK(41,21) +72(41,21) 
DIMENSION S$(414521)2+VOR(41.21) 5sFJ5(41) 

DIMENSION X(41)+RA1(°41)+.SWT( 200) 

DIMENSION UC41521)0V (41521 ),W(41521) 

DIMENSION GZW(41+.21)eWT(41.21) »SDZ(200) 

DIMENSION XT(41 521) 

COMMON /C1/ ZZ(1500)sDXsKZ > BI 
COMMON /C2/ A(41s21)2BsC(41521)50(41521)sEsF( 41521) 


1,WW(41521) 


COMMON /C4/ Mi1.N1 

COMMON /C3/ MoN 

COMMON /C6/ AREA 

REAL MAXS»MAXT 

READ(5510) MeaNoR 

WRITE(6210) MsNeR 
FORMAT(2155F5e1) 

READ(5s214) PR»sRC»REYsDEANSBI 
FORMAT(5F15.5) 

WRITE(6514) PRsRCsDEANsREY 
WRITE(65114) BI 
FORMAT(*BIGT NUMBER =!,F25.3) 
FORMAT (//s"*PRANDTL NUMBER ="5F15.59/5 


1 "CURVATURE ='5F15<.5s/% 
* *DEAN NUMBER =* 6F15.59/7s*REYNOLD NUMBER ="'5F15.5) 


M1=M+1 

Ni=N+1 

MM 1=M~—1 

NM1=N-1 

HI=(1.t+1./R)/(2.*M) 

HJ=(1 o+R)S(4.N) 

DX=HI/2. 

DX 2=HI*HI 

DO 450 I=1.M1 

READ(5s460) (SC IsJ) sJ=19N ) 
DO 451 I=1.M!1 

READ(5s460) (WCIsJ) sJ=15N ) 
DG 452 I=1.Mi1 

READ(53460) (VOR(IsJ)sJ=15N1) 
DO 888 I=1.Ml1 

READ(5s460) (UC IsJ)sJ=1.N1) 
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DO 889 TI=1.Mi1 
READ(52460) (VC IoJ) sJ=19N1 ) 
FORMAT(5E16.7) 
ENTRANCE CONDITICNS 
DO 15 I=1,.M1 
w(I»«N1i)=0. 
DG 15 J=i1,N1 
FT Ol ss =1 76 
DO 100 I=1.M1 
X (1 )=-(1.41./R)44.4( 1-1) *HI 
RAIC I )=1.4+X(1)/RC 
DO 12 I=1>5Ml1 
DO FP2/J=11N1 
WC IT sJ)=W(I oJ)/SREY 
DO 18 I=1.M1 
DO 18 J=1,N1 
WW(C To J) =2.*W( I sd) *DX2/RAI (I) SPR 
PR1=1./PR 
DO 11 I=1.M1 


dO 11 J=1sN1 


ACI sJ)=PR1I-DX* (UC Ts J)—-1./ 0 RC#XRAI(I)) 
B=-2.*PR1 
C(IsJ)=PR14DX*¥(UC Ts J)—-1-/°RC*RAI(I)) 
OC I.J)=—-PR1+DX*V(IoJ) 
E=2.*PR1 
F(IsJ)=-PRI-DX*V(IsJ) 
CONTINUE 
AREA=(1.4R)*(1.4+1./R)/36 
KMAX=1000 

Z=0% 

DX=HI 

Z=0% 
KZ=0 
KK=1 

KZ=KZ+1 

DO 19 I=1.M1 

DO! £O Jalen 
TK(CIsJ)=T(IsJ) 

AXIAL STEP 

DZ=5.E-5 

IF (KZ.GT~100) OZ=1.E-4 
IF (KZ.eGT.450) DZ=2.E-4 
IF (KZ.GT-650) DZ=S5.E-4 
IF (KZ.eGT.850) DZ=1.E-3 
Z=Z+DZ 

ZZ(KZ)=DZ 
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WRITE(62201) KZsZ 

DO 22 I=1>5Mi 

HGY22"3=1.N1 

TK (CIT eJS)=TCI5J5) 

CALL TEMP (T.TK) 

DO 23 I=1.M1 

DO 23 J=1,N1 

TICIsJ)J=T(CI oJ) *¥WCI 2S) 

CALL SIMPS(T1sMsNsHIsHJsSUMM) 

TBK=SUMM 

SUT=0. 

DO 42 I=1.M1 

SUT=SUT4#T(I5N1) 

DO 44 J=1.5N 

SUT=SUT#T(1sJ)4+T(M15J) 

SUT=2.*SUT-T(151)-T(M1 91) 

AVET =SUT/(2*(M4+2*N) ) 

TAW=AVET 
ANU1=-BI*TAW/(TAW-TBK) 


CALCULATIONS OF TEMP. GRADIEND IN FLOW DIRECTION 


DO 25 J=1>.N1 

DO 25 I=1.Mi1 

GZW( Ts JJ =H=CTOCI J) —-TKC Id) ) *WOIs I) /ZZ(KZ) 
CONTINUE 

CALL SIMPS(GZWseMsNsHIsHJs SGWT) 

ANU2 = C.25% SGWT/ (TAW-TBK) 
AVNU=0 .5* (ANU1+ANU2) 

DNU=100.* ABS(ANU1—ANU2)/AVNU *0.5 

WRITE(6,501) ANUI1 » ANU2 s AVNUsDNU 
Tw=T(1.N1) 

SUPER=(1.-TW)/TW 

WRITE(6,502) TSK »sTAW ».SUPER 
KKK=KZ/30 

IF (KKK.NE.KK) GO TO 500 

KK=KK+41 

WRITE(62503) 

DO 161°E8=145 M1 

WRITE(63504) I 

WRITE(65505) (TCIsJ) sJ=1.Ni1 ) 

FORMAT(/s*NO. GF SECTION®+.15510Xs*DISTANCE FROM 


1 ENTRANCE® .F10.6) 


FORMAT (30X » *¥ **",3F15.62E20.5) 
FORMAT ( *NU1%,F10.6910Xs *NU2Z" 5F1066910Xs "AVERAGE 


1NU *sF10.6310X3E20.4) 


FORMAT(*®*BULK TEMP." s £15.5 s*WALL TEMP.%sE15.5s 


1 "SUPERHEATED RATIO .E15.3) 
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FORMAT(// *TEMP. DISTRIBUTIONS*/) 
FORMAT (*I* ,15) 

FORMAT( © €,9£12.4) 

CONTINUE 

IF (KZ.GEe«KMAX) GO TO 90 

GO TO 20 

CONTINUE 

STOP 

END 

SUBROUTINE TEMP(T TK) 


‘DIMENSION 1041221) sTK( 41921 )s TKH( 41521) 5X(41) 


rey 


DIMENSION A1(41)sA2(41)+A3(41)5A4(41)5A5(41) 
COMMON /C1/ ZZ(1500) sDXsKZ » BI 

COMMON /C2/ A(41.21).B.C(41.21),0(41521) sEsF( 41521) 
oWW(41,421) | 

COMMON /C3/ MsN 

COMMON /C4/ MiNi 

DO 20 J=15N 

DO 10 I=25M 

A1(1)=CCI.sJ) 

A2(1)=B-wwCI »J)/ZZ( KZ) 

AZ(I)=ACI.J) 

IF (J-EQ.1) GO TO 15 

AG (IT)=D(I1oJS)*TKC Io S41) +CE-WWOI oJ) ZZ0KZ)) *TKC ISS) + 


1F(IsJ)*TK(CIsJ5-1) 


GO TO 10 
A4(I)=COCI eS) FEC Io Jd) XTKC Io Jt1)4+(CE-WW( 1 oJ) /ZZ0KZ) )* 


1TK(IsJ) 


CONTINUE 
A2(1)=-2. -2.*DX*BI 
A3(1)=2. 
A1(M1)=2. 
A2Z(M1 )=-2.-2-*DX*BI 
IF(J.EQ.1) GO TO 16 
AS (1) =—-TK(C 1.541 )426*TK(C 155) -TK(15/d-1) 
AG(M1 )=—-TK(M1 > J41)420*TKOCM1 95) -TKOM1,J—-1) 
A4( 1) =-2.*TK(1 52)42-*TK(1 91) 
A&4 (M1) =—-2e*TK(M152)4+2.*TK( M151) 
CALL GAUSS(A1,A2sA3sA45M1 5X) 
DO 17 I=1.M1 
TKH(IsJ)=X(I) 
CONTINUE 
DO 21 I=1>,Ml1 
A1(1I)=1. 
A2(1)=-4. —-2.*DX*BI 
A3(I)=1. 
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A4(I )=-2.*TKHCI.N) 

A1l(M1)=2. 
A2(1)=A2€1)-2.*DX*BI 

A2 (M1 )=A2€M1)-2.*DOX*BI 

CALL GAUSS (AilsA2sA3sA49M1 9 X) 
DO\22 (1=13M1 

TKHCI5N1)=X(1) 

DG 40 I=25M 

DO 30 J=1.N 

AL(JJ=FC 155) 
A2(J)=EtWWCI J) /ZZ(KZ) 
A3(J)=D0(C1.J) 

AG (J)=AC Is JI) *¥TKHCI41 55) 4° BtWWw( I sJ)/ZZ0KZ) ) * 


1TKHCIsJ)4+CCIsJS) *TKHC I-14) 


CONTINUE 

A1(N1)=2. 
AZ2(N1)=-2. —-2.*DX*BI 
A3(1)=F(151)4D( 151) 
A4 (Ni )=—-TKH(CI419N1)+2.*TKHCI9Ni )-TKH(I-15N1) 
CALL GAUSS(A1sA29A3 sA45Ni1 5 X) 
DO 40 J=1,N1 
TCLs SD=XOS) 
DO 50 J=1.N1 

A1l(J)=1. 

A2(J)=-4. —-2.*DX*BI 
A3(J)=1. 

A4 (059) =-226*T( 255) 

AS (J) =—-20*T(M5 J) 

A3(1)=2. 

A1(N1)=2. 

A2(N1)=A20N1 )-2.*DX*BI 
CALL GAUSS (A1lsA2sA3sA4sN15 X) 
DO 52 J=1,5N1 
TC1sJ)=X(J) 

CALL GAUSS(A12A2sA3 sASsN1 9X) 
DO 54 J=1.N1 
TOM1sJ)=XCJ) 

RETURN 

END 


SUBROUTINE FOR GAUSSIAN ELIMINATION 
SUBROUTINE GAUSS(As85CsDsN1 9X) 
DIMENSION A(61)+8(61)5C(61) »0(61)+X(41).BETA(61) >» 
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1GAMMA(61) 


N=Ni-1 
BETA(1)=B(1) 


- GAMMA(1)=D(01)78(1 ) 


DO FPO Lt=2 sN1 

BETACI )=B(1)-AC 1) *CCI-1)/BETACI-1) 
GAMMA (I )=(D(1I)—A(C 1) *GAMMAC I-1))/BETA(I1) 
X(N1)=GAMMACN1 ) 

T=N1 

DO 20 K=1.5N 

I=I-1 

XC I )=GAMMACI)-CC IT) *XCI41)/BETACI) 
RETURN 

END 


SIMPSON"S RULE FOR SINGLE INTEGRATIGN--(FCJ) »J=159N1) 
SUBROUTINE SIMP(FsNsHJ»sSUM) 
DIMENSION F(61) 

Ni=N+1 _ 

NM1=N-1 

SUM=F (1)4+FCON1) 

DO 10 J=2seNo2 
SUM=SUM+4 4. 0*F( J) 

DO 12 J=3sNM15.2 
SUM=SUM+#+2.0 *F(J) 
SUM=SUM*HJ/3. 

RETURN 

END 


SIMPSON"S RULE FOR DOUBLE INTEGRATION--CCFC(IsJ)5 


1J=1.oN1) »T=1.M1) 


SUBROUTINE SIMPS(FsMsNsHIsHJs SUMM) 
DIMENSION SUM(61) .F (41.41) 

Mi=M+1 

N1=N+1 

MM 1=M—1 

NM1=N-1 

DO 8 I=1>5M1 

SUMCI)=FC1I51)4F(1I2N1) 

DO 4 J=29No2 

SUM(C I )=SUM(1)4+4.04F(I5J) 
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DO 6 J=3eNM1.2 
SUM(I)=SUM(1I)+2.0%F (14 J) 


SUM(I)=HJ*SUM(1I)73. 


SUMM=SUM(1)+SUM(M1 ) 
DO 10 I=2sMs2 
SUMM=SUMM4#4,.0*SUM(I1) 
DO 12 I=35MM1.2 
SUMM=SUMM+42.0*SUM(CTI ) 
SUMM=HI*SUMM/3. 
RETURN 

END 
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NUMERICAL RESULTS IN THERMAL ENTRANCE REGION WITH 


UNIFORM WALL TEMPERATURE 


= 0.7 , K = 368 
Zz Nu, 

5x10°> 14.748 
60 13.694 
70 12.953 
1x107¢ 11.576 
12 10.983 
15 10.323 
20 9.578 
25 9.092 
30 8.765 
35 8.544 
40 8.399 
50 8.261 
60 8.258 
70 8.339 
80 8.473 
90x10" 8.641 
1x10"! 8.828 
12 9.230 


APPENDIX C.1 
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NUMERICAL RESULTS IN THERMAL ENTRANCE WITH 
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